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AN IMPROVED VERSION OF THE NASA-LOCKHEED 
MULTIELEMENT AIRFOIL ANALYSIS COMPUTER PROGRAM 

G. W. Brune and J. W, Manke 
Boeing Commercial Airplane Company 


SUMMARY 


^ This document contains a description of an improved version of the NASA-Lockheed 
computer program for the performance prediction of high lift airfoils. Modifications of 
the aerodynamic model and the computer program include: 

• Boundary layer and wake displacement effects are represented by an equivalent 
distribution of sources along the airfoil surface and along the wake centerlines. 

• Wake parameters are predicted using the lag-entrainment method of Green. 

• Profile drag is calculated by the Squire and Young formula. 

• Para mete rs of ordinary turbulent boundary layers are calculated by the method of 

Nash and H icks. [ _ _ 


• Onset of the separation of confluent turbulent boundary layers is determined by a 
modification of Goradia’s confluent boundary layer method. 

• High lift airfoils with up to 10 components can be analyzed. 

• The Boeing version of the computer code is well structured featuring new control 
routines and new subroutines for geometry and potential flow calculations. Old 
subroutines are thoroughly commented and documented. 

The program is evaluated by comparison with recent experimental high lift data, 
comprising lift, pitching moment, and profile drag, as well as, detailed distributions of 
surface pressures, boundary layer integral parameters, skin friction coefficients, and 
velocity profiles. The results of this evaluation show that the contract objectives of 
improving program reliability and accuracy have been met. 



INTRODUCTION 


HISTORICAL REMARKS 

In the past, high lift design and technology rested in the hands of a few experienced 
aerodynamicists. Design methodology and criteria were heavily influenced by the 
analytical inviscid flow methods and the experimental data available. With the advent 
of high-speed computers and the appearance of improved models for turbulent flows, 
many of these complex problems, including high lift design and analysis, were attacked 
theoretically. 

One such approach to high lift or multielement airfoil analysis was developed by 
Goradia and his coworkers (ref. 1) at Lockheed-Georgia under the sponsorship of the 
NASA-Langley Research Center. This program was among the first attempts at 
analyzing the complex viscous flow about slotted airfoils and has received worldwide 
distribution and usage. A unique feature of this multielement airfoil program is the 
model of the confluent boundary layer flow (ref. 2). 

Over the years, the original version of the program was modified extensively to improve 
its predictions for different types of high lift airfoils. Many improvements, mainly in the 
area of the potential flow calculation, were made by researchers at the Langley 
Research Center (ref. 3). For this reason, the code is generally referred to as the 
NASA-Lockheed multielement airfoil program. A version for single element airfoils was 
recently extracted from the multielement airfoil code by researchers at North Carolina 
State University (ref. 4). 

Widespread and steady usage of the computer program clarified its strengths and 
weaknesses. Both favorable and unfavorable aspects have been brought to the surface 
by continued attempts at using the program as an engineering tool. The more serious 
shortcomings wcto the lack of agreement between the documentation and the available 
version of the code and the high-failure rate in applying the method for various 
configurations. However, the program was found to contain sufficient positive features 
to justify its choice as a starting point for future theoretical work in the high lift area. 

In July 1976 work was begun in the high lift research group of The Boeing Company on 
a joint program with NASA-Langley to evaluate and improve the NASA-Lockheed 
multielement airfoil code. The work consisted of two phases. The first phase had the 
objective to document and evaluate the "baseline” version of the code which was 
supplied to Boeing by NASA-Langley prior to the beginning of the contract work. In 
addition, certain minor improvements of the aerodynamic methodology were made 
during the first work phase. In February 1977, the phase one version of the computer 
code was delivered to the Langley Research Center together with a detailed 
documentation of its underlying aerodynamic theory (ref. 5). 



The second phase of the contract work involved a major revision of the flow model used 
in the NASA-Lockheed program that in turn required substantial modifications of the 
computer code. This document contains the results of the second phase of the contract 
work including the complementary Boeing IR&D work. The evaluation of the 
predictions of the various versions of the computer code by comparison with recent 
experimental data of high lift airfoils is described in a separate document (ref. 6). 
However, this document contains a few of these comparisons to provide the user with a 
reasonably self-contained guide to the latest version of the computer program. 

MULTIELEMENT AIRFOILS 

The flow around high lift airfoils is characterized by many different inviscid and viscous 
flow regions. Their complex physics is illustrated in figure 1. In particular, the 
existence of confluent boundary layers and the regions of separated flow distinguish 
the high lift airfoil problem from the aerodynamic problem of airfoils at cruise 
conditions. The various flow regions, including the outer potential flow, the ordinary 
laminar and turbulent boundary layers, viscous wakes, and the confluent boundary 
layer, are analyzed by the code. Furthermore, the prediction of transition from laminar 
to turbulent boundary layer flow and the prediction of the onset of boundary layer 
separation are a necessary part of the code. Cove separation and large scale separation 
phenomena, however, are not modeled. 

The computer code described in this document calculates the flow about high lift airfoils 
assuming that: 

• The flow is attached to the airfoil's surface 

• The flow is two-dimensional and subcritical 

• The high lift airfoil consists of up to 10 components 

These are the main assumptions of the multielement airfoil method. Additional 
assumptions are listed in the pertinent sections of this document. 

MODIFICATIONS OF THE AERODYNAMIC MODEL 

The described aerodynamic theory differs from the theory of the baseline version of the 
NASA-Lockheed program in the following areas: 

• The method used to represent the effect of the viscous flow on the outer potential 
flow within the overall iterative solution procedure, that is termed equivalent 
airfoil representation in the baseline version of the program, has been modified. It 
has been replaced by the surface transpiration method which uses an equivalent 
distribution of sources along the airfoil surface and the wake centerlines to model 
the boundary layer and wake displacement effects. 
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• The flow model of the potential core region has been changed. The new method 
performs independent boundary layer and wake calculations. These calculations 
utilize the ordinary laminar and turbulent boundary layer routines of the baseline 
version of the code, and in addition, the lag-entrainment method of Green (ref. 7) 
for wake flows. The revised flow model of the core region calculates the location of 
the wake centerlines. 

• An attempt is made to predict the onset of separation of the confluent boundary 
layer by a modified version of Goradia’s confluent boundary layer method. In this 
method, the power law velocity profile of the wall layer is replaced by Coles' 
two-parameter velocity profile (ref. 8). The latter is known to be an adequate 
representation of thin turbulent boundary layers. 

• The drag prediction method of Squire and Young (ref. 9) has been incorporated into 
the program, replacing the previous pressure and skin friction integration scheme. 

• The original method used for the prediction of separation for ordinary turbulent 
boundary layer flow has been replaced by the Boeing version of the method of Nash 
and Hicks (ref. 10). 

• The slot flow calculation has been removed from the program. The flow in the slot 
between adjacent airfoil components is now calculated as an integral part of the 
overall computation. This is based on the same potential flow and viscous flow 
models, including the same account of compressibility effects, that are used in the 
remainder of the flow field. 

• The code has been made operational for negative airfoil overlap. 

• Several logical errors in the formulation of the aerodynamic model and its 
numerical implementation have been corrected, 

I COMPUTER CODE | 


The outlined modifications of the aerodynamic theory required a major overhaul of the 
computer code. Most parts of the code have been rewritten using a systematic approach 
to computer software design. This work was guided by a functional decomposition of the 
many aspects of the aerodynamic model and its numerical implementation. In addition, 
a detailed study was made of the data flow within the program, and the logic of the code 
was outlined prior to the actual program development using a pseudo code. The most 
important results of this work, such as the higher levels of the functional 
decomposition, a brief description of the data structure, and a unified list of symbols, 
are included in this document. 

All control routines, the geometry package, and the potential flow routines of the 
program have been replaced. Other subroutines performing such functions as tracing 
streamlines, computing wake flow characteristics, predicting confluent boundary layer 
separation, etc., have been added to the code. However, several major subroutines of the 



old code including LAMNA, TURBL, TURB, C0NF7, CONF8, and some mathematical 
subroutines, had to be retained with only minor modifications. The schedule of the 
contract did not permit time to restructure these routines. Nevertheless, a considerable 
amount of time was spent reorganizing the COMMON blocks of these old subroutines, 
rationalizing conflicts of symbols, and interfacing the data structure with the new code, 
A list of symbols and comment cards, sufficient to guide an experienced user through 
the code, were added to each of the old subroutines. 

The new routines in the computer program use dynamic data structures that do not 
limit explicitly the number of surface points representing the airfoil geometry. The 
input format and the boundary layer routines of the baseline version of the code limited 
the number of computational surface points to 165 and the number of data points to 65 
per upper or lower airfoil surface. 

Previous users of the computer program should note that small changes have been made 
in the input and output formats. 


ON THIS DOCUMENT 

The document combines the results of contract work and complementary Boeing IR&D 
work. The Boeing IR&D funds supported the following work and documentation. 

• The confluent boundary layer model. 

• The top-down design of the computer code including the functional decomposition 
of the aerodynamic theory, the investigation of the data flow, and the pseudo code. 

• The evaluation of the computer program by comparison with experimental data of 
McGhee and Beasley; Wentz, Seetharam, and Fiscko; Foster; Ljungstrbm; and The 
Boeing Company. Most of these results are described in a separate document 
(ref. 6). 

The table of contents closely follows the functional decomposition of the code. The 
letters behind a heading refer to the corresponding function of the functional 
decomposition. No attempt has been made to document the code on a subroutine by a 
subroutine basis. 

Those sections of the document describing the theory of old subroutines have an 
individual list of symbols. All other sections share one common list of symbols relating 
the theory to the computer code. Cross referencing from section to section has been 
avoided. 
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STRUCTURE OF THE COMPUTER CODE 


The described version of the NASA-Lockheed multielement airfoil computer program 
conforms to the Langley Research Center computer programming standards. It is 
Avritten in the CDC FORTRAN Extended 4 (FTN4) language and will run under the 
CDC Network Operating System (NOS). Program I/O is performed only by FTN4 
statements using the standard system file names INPUT (TAPE5) for card reading and 
OUTPUT (TAPE6) for printing. 

In the sections below, the design of the computer code is discussed, the overlay structure 
of the code is described, and a short description of each subroutine is given. 

DESIGN OF THE CODE 

The programming methodologies used to design and develop the new version of the 
computer code include: 

• Functional decomposition 

• Data flow analysis 

• Control flow analysis 

Each of these interrelated design tasks were performed several times in an iterative 
manner to produce a final design for the new version of the computer code before any 
changes or improvements to the baseline code were made. The final design for the new 
version of the code resulted in major changes in the three following program sections; 
upper level control routines, geometry preprocessing routines, and the potential flow 
solution routines. The final design was also used to integrate the major new 
aerodynamic models into the baseline code; they include the representation of the 
displacement thickness with sources, Green’s wake solution technique, and the modified 
confluent boundary layer method. Table 1 lists all subroutines in the baseline version of 
the code and indicates the type changes made to incorporate them in the new version. 

The functional decomposition of the code was based on the engineering specification of 
the aerodynamic models and the numerical techniques necessary for their solution. 
Figures 2 and 3 show the upper level decomposition charts where the major functions 
are defined in engineering terms. The complex physics of the flow about multielement 
airfoils is reflected in these charts, e.g., the potential flow solution, the ordinary 
laminar and turbulent boundary layer solutions, the confluent boundary layer solution, 
and the wake layer solution. The charts also demonstrate that the overall iteration of a 
solution is a high level function; and, the prediction of the transition from laminar to 
turbulent boundary layer flow and the prediction of the occurence of separation are 
functions of importance equal to the laminar and turbulent flow solutions. All new 
subroutines in the code identify which module in the functional decomposition they 
implement. 



Table 1. — Modifications of the Subroutines in the Baseline Version 


Subroutine 

No 

change 

Minor 

change 

Major 

change 

Delete 

MAIN 



X 


POINT 

X 




SLOPE 

X 




TRANS 



X 


DISTP 

X 




FTLUD 




X 

DIR 

X 




LSQ 




X 

PROOT 

X 




MA1N1 



X 


READIT 


X 



GEOM 



X 


ROTRAN 



X 


AS LOT 



X 


NORMAL 



X 


MAIN 2 



X 


CHEN 



X 


MATRIX 



X 


POTLF 



X 


CAMBER 




X 

SMOOTH 

X 




VOVBT 




X 

THICK 




X 

COMPR 



X 


STAG 



X 


MAIN 3 



X 


LOAD 



X 


LAMNA 


X 



BLTRAN 


X 



TURBL 


X 



TURB 


X 



DERIV 


X 



START 


X 



CONFBL 




X 

CONF 5 




X 

CONF 7 


X 



CONF 8 


X 



DLIM 

X 
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Figure 2. — Functional Decomposition of NASA- Lockheed Program 





Figure 3. — Functional Decomposition of Viscous Flow Solution 

















The data flow analysis of the code was done with the aid of HIPO charts: an example of 
a completed HIPO chart is given in figure 4. For each module identified in the 
functional decomposition, the input and output data for the module, as well as its 
decomposition, are specified on the chart. Control of the data flow within a module is 
maintained by requiring that the input to any of its submodules must be either an input 
to the module or the output of another of its submodules. Once HIPO charts have been 
completed for all modules in the functional decomposition, all data groups have been 
identified, and the data flow specified. The data groups identified for the new version of 
the code by this process are listed in the section of this document titled Symbols. 

The control flow analysis of the new code was done with the aid of pseudo code; an 
example is given in figure 5. Pseudo code is a small set of simple logic and loop 
statements which suffice to describe the control within a module of the functional 
decomposition. Although the submodules of a particular module can be used in any 
sequence and any number of times to complete the function of the module, it is an aim 
of the design process to keep the control within a module as simple as possible. All new 
subroutines in the code include as comment cards the pseudo code for the module which 
they implement. 


OVERLAY STRUCTURE 


The CDC overlay system is used to assure that the computer code will execute in a field 
length less than 100 K (octal). The division of the code into overlay sections follows the 
decomposition of the solution process: user input processing and geometry data 
preprocessing, the potential flow solution, and the viscous flow solution. The viscous 
flow solution divides into the laminar flow solution, turbulent flow solution, wake layer 
solution, the confluent boundary layer solution of Goradia, and the modified confluent 
boundary layer method. The overlay structure of the program is described in detail by 
the table in figure 6. 

I SUBROUTINE description! 

This section contains a description of those subroutines which perform the analysis. The 
subroutines can be divided into several groups according to their function; major control 
routines, user input and geometry data processing, potential flow solution, viscous flow 
solution, and library routines. In the succeeding descriptions the subroutines are 
divided into these groups. 

The execution of the program is directed by the major control routines which are listed 
below: 

MAIN 

SYSOL 


INITR 


Provides the primary logic and data flow control for the entire program. 

Provides the logic and data flow control for the solution iteration for each 
requested angle of attack and freestream Mach number. 

Initializes parameters of the solution iteration procedure. 
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Figure 4. — Sample HlPO Chart 




BEGIN MODULE BBd 

FOR EACH AIRFOIL COMPONENT DO 

FOR EACH AIRFOIL COMPONENT SURFACE W 

BEGIN MODULE BBCA : COMPUTE LAMINAR BOUNDARY 

LAYER solution 

If. TRANSITION TO TURBULENT FLOW THEN 

BEGIN MODULE BBCB : COMPUTE TURBULENT 

BOUNDARY LAYER SOLUTION 

ENDIF 

ENDDO 

ENDDO 

FOR EACH AIRFOIL COMPONENT M 

BEGIN MODULE BBCC : COMPUTE WAKE LAYER SOLUTION 

]F NOT THE LAST AIRFOIL COMPONENT AND CORE REGION 
ENDS BEFORE THE TRAILING EDGE THEN 

BEGIN MODULE BBCD : COMPUTE CONFLUENT BOUNDARY 

LAYER SOLUTION FOR NEXT AFT 
COMPONENT 

I£ CONFLUENT BOUNDARY LAYER FLOW SOLUTION 

DEGENERATES INTO ORDINARY TURBULENT FLOW THEN 

BEGIN MODULE BBCB : COMPUTE TURBULENT BOUNDARY 

LAYER SOLUTION 

ENDIF 

ENDIF 

ENDDO 

END MODULE BBC 


Figure 5. — Sample Pseudo Code 



OVERLAY (0,0) 





PROGRAM : MAIN 





SUBROUTINES : SYSOL 

LOAD 

srfit 

REWOX 

PROOT 

INITR 

DYNSET 

SMOOTH 

SUMRX 

GLESOS 

SOLVR 

GETBDA 

initx 

REAOR 

FBSUBS 

CONVR 

PRTBDA 

READX 

blkit 

OECOM 

SRCER 

FREBDA 

writx 

ERSET 

VIPD 




OVERLAY (1.0) 

PROGRAM : HAIN10 

SUBROUTINES : 

INPTR LOFTR 

READIT SLOTR 

GEOH GLOBD 

GEOMA GEOMC 

RESCL ANGLR 

AFPRM LOCLD 

COMPT SURFD 

SMOPT TRANF 

GEOMB DISTP 


DIR 

WAKCL 

FTLUD 


OVERLAY (3,0) 
PROGRAM : MAIN 30 


SUBROUTINES 

POTLF V 
POTLFA V 
POTLFB F 
POTLFC C 
POTLFO 5 
POTLFE / 
ANLYS ( 
WAKET ; 
WAKEJ ( 


WAKES 

WAKEG 

NEWTR 

COMPR 

STAGN 

AlCVAL 

CPTAIC 

SAVA I C 

GETAIC 


OVERLAY (4.0) 

PROGRAM : MAIN40 

SUBROUTINES 
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SOLVE 


Provides the logic and data flow control for one step in the solution 
iteration procedure: this includes representation of the boundary layer, 
potential flow solution, viscous flow solution, and loads estimate. 

CONVR Checks the convergence criteria to determine when the solution iteration 

procedure can be terminated. 


Listed are the subroutines that read and analyze the user input data and preprocess the 
geometry. 


MAINIO 

INPTR 

READIT 

GEOM 

GEOMA 

RESCL 

AFPRM 

COMPT 

SMOPT 

GEOMB 

LOFTR 

SLOTR 


Provides the primary logic and data flow control for reading the user 
input and preprocessing the geometry. 

Provides the interface between the data structures of the new version of 
the program and READIT, the input reading routine of the baseline 
version. 

Reads the user input data cards and stores the problem description. Some 
format and consistency checking of the input is performed. 

Provides the primary logic and data flow control for the preprocessing of 
the geometry data. 

Provides the logic and data flow control for the airfoil parameter 
determination phase of the geometry analysis. 

Stores the input geometry data in the internal basic data array format, 
and applies the user specified geometry scaling factor. 

Determines the basic airfoil parameters; e.g., number of computational 
surface points, airfoil chord length. 

Computes the computational surface points. The input surface points are 
redistributed by an algorithm based on curvature (DISTP). 

Smooths the geometry data computed by COMPT with a simple 
smoothing algorithm (SMOOTH). 

Provides the logic and data flow control for the lofting of the geometry 
data in the global coordinate system. 

Components other than the main component, are rotated and translated 
into the coordinate system of the main component. 

Analyzes the slot geometry of the lofted airfoil. 
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GLOBD Computes the initial wake centerline geometry and stores the final lofted 

airfoil geometry data in the global coordinate data array. 

GEOMC Provides the logic and data flow control for the calculation of the local 
coordinate geometry data. 

ANGLR Calculates the global-to-local coordinate system transformation data. 

LOCLD All components are rotated and translated into their respective local 

coordinate system. 

SURFD Computes the surface fitted coordinates (arclength) for each component. 

WAKECL Computes the initial wake centerlines. 

The subroutines that calculate and analyze the potential flow solution are listed below: 

MAIN30 Provides the primary logic and data flow control for the calculation and 
analysis of the inviscid flow solution. 

POTLF Provides logic and data flow control for calculation of the incompressible 

surface velocity by Oeller’s potential flow solution technique. 

POTLFA Calculates the solution matrix determined by Oeller’s method and the 
specification of the Kutta condition. 

POTLFB Provides an interface with a linear equation solution package and checks 
the numerical conditioning of the solution matrix. 

POTLFC Calculates the right hand side determined by the freestream velocity, 
wake centerline sources, and the Kutta condition. 

POTLFD Provides an interface with a linear equation solution package and 
computes the vortex strengths determined by the right hand side. 

POTLFE Calculates the incompressible surface velocity from the freestream 
velocity and the source and vortex strengths. 

ANLYS Provides the logic and data flow control for the analysis and correction of 

the incompressible surface velocity. 

WAKET Computes the initial values of the wake centerline parameters for the 
update of the wake centerline geometry. 

WAKEJ Computes the Jacobian matrix for the wake centerline parameters for the 

update of the wake centerline geometry. 

WAKES Computes the stream function value of the potential flow solution 

velocity field at the corner points of the wake centerline. 
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WAKEG Computes the geometry of the updated wake centerline from the values 
of the wake centerline parameters. 

COMPR Applies the compressibility corrections to the incompressible surface 

velocity and computes the local Mach number and pressure coefficient. 

STAGN Estimates the location of the stagnation point of the flow field for each 

airfoil component. 

The subroutines which calculate the boundary layer parameters for the viscous flow 

solution are: 

MAIN40 Provides the primary logic and data flow control for the calculation and 
analysis of the viscous flow solution. 

VSFINA Provides an input interface between the data structures of the new 
version of the program and the laminar and turbulent boundary layer 
analysis routines of the baseline version. 

VSFINB Provides an input interface between the data structures of the new 
version of the program and the confluent boundary layer analysis 
routines of the baseline version. 

LAMNA Computes compressible laminar, boundary layer flow solution using an 

integral method similar to the method of Cohen and Reshotko. 

BLTRAN Computes transition of compressible laminar flow to turbulent flow or 
separation of the laminar boundary layer. 

TURBL Computes incompressible, turbulent, boundary layer flow solution using 

an integral method similar to the method of Truckenbrodt. 

TURB Computes incompressible, turbulent boundary layer flow solution using 

an integral method due to Nash and Hicks. The method can predict 
separation of the turbulent boundary layer. 

DERIV Computes values of partial derivatives of parameters defined in the Nash 

and Hicks method. 

START Computes initial values for the parameters defined in the Nash and 

Hicks method. 

WAKEI Computes the initial values of the parameters for Green’s wake layer 

solution method. 

WAKED Computes derivatives of the parameters for Green’s wake layer solution 
method. 
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WAKEP 

ECORE 

CONF7 

CONF8 

DLIM 

CONF 

CONFIl 

CONFDl 

CONFPl 

CONFI2 

CONFD2 

CONFP2 

VSFOUA 

VSFOUB 

SRCER 


Computes the values of the parameters for Green’s wake layer solution 
method at each wake centerline corner point. 

Estimates the end of the potential flow core region in the flow field 
behind each slot. 

Computes and displays the parameters of Main Region I in the analysis 
of the confluent boundary layer with the method due to Goradia. 

Computes and displays the parameters of Main Region II in the analysis 
of the confluent boundary layer with the method due to Goradia. 

Limits the magnitude of estimates of derivatives used in CONF7, 
CONF8. 

Provides the logic and data flow control for the analysis of the modified 
confluent boundary layer method. 

Computes the initial values of the parameters of Main Region I in the 
analysis of the modified confluent boundary layer method. 

Computes the derivatives of the parameters of Main Region I in the 
analysis of the modified confluent boundary layer method. 

Computes and displays the parameters of Main Region I in the analysis 
of the modified confluent boundary layer method. 

Computes the initial values of the parameters of Main Region II in the 
analysis of the modified confluent boundary layer method. 

Computes the derivatives of the parameters of Main Region II in the 
analysis of the modified confluent boundary layer method. 

Computes and displays the parameters of Main Region II in the analysis 
of the modified confluent boundary layer method. 

Provides an output interface between the data structures of the new 
version of the program and the laminar and turbulent boundary layer 
analysis routines of the baseline version. 

Provides an output interface between the data structures of the new 
version of the program and the confluent boundary layer analysis 
routines of the baseline version. 

Computes the source strength values to represent the displacement 
thickness of the viscous layers. 
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LOAD 


Computes the global aerodynamic parameters: lift coefficient, drag 
coefficient, pitching moment coefficient, axial - and normal force 
coefficients. 

The library routines support the control and analysis routines listed above. Among the 
functions supported by the library routines are: dynamic storage control, aerodynamic 
influence coefficients (AIC’s) calculation, solution of a system of linear equations, 
Newton algorithm, and integration of a system of ordinary differential equations. The 
library routines are: 

DYNSET Initializes and controls the dynamic storage work area. 

GETBDA Reserves storage in the dynamic storage work area for one of the several 
types of basic data arrays. 

FREBDA Frees the storage in the dynamic storage work area which was assigned 
to a basic data array by GETBDA. 

PRTBDA Displays the contents of a basic data array which was assigned by 
GETBDA. 

AICVAL Provides the logic and data flow control for the calculation of the 
aerodynamic influence coefficients (AIC’s) for Oeller’s potential flow 
solution technique. 

CPTAIC Computes the stream function and velocity AIC’s for vortex and source 

distributions on the surface (and wake) segments. 

SAVAIC Saves the AIC’s generated by CPTAIC on an I/O unit. 

GETAIC Reads the AIC’s stored on an I/O unit by SAVAIC. 

GAUSSR Provides an in core, Gaussian, linear equation solution package with a 
pivoting capability. 

RSEITE Provides a second solution capability for the algorithm implemented in 

GAUSSR. 

REDUCE Provides an out of core, Gaussian, linear equation solution package with 
a pivoting capability which reduces the system to triangular form. An in 
core equation solver is available upon request. 

BCKSUB Completes the out of core Gaussian solution process in REDUCE by doing 
the back solution process. 

NEWTR Provides the logic and data flow control for solving a system of nonlinear 

equations with a Newton algorithm. 
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INTRG 

TRANF 

DISTP 

DIR 

POINT 

SLOPE 

PROOT 

SRFIT 

SMOOTH 

INITX 

READX 

WRITX 

REWDX 

SUMRX 

READR 

BLKIT 


Provides the logic and data flow control for integrating a system of 
ordinary differential equations with an integration algorithm. 

Rotates and translates a coordinate geometry array. 


Computes a coordinate geometry array such that the corner points are 
separated by equal increments of curvature of the surface. 

Computes an estimate of the derivative of a tabulated function at one of 
the tabular points. 

Computes an estimate of the value of a tabulated function at a value of 
the independent variable with parabolic interpolation. 

Computes an estimate of the derivative of a tabulated function at a value 
of the independent variable with parabolic interpolation. 

Computes the roots of a cubic polynomial. 

Uses a parabolic fit to surface coordinate points to estimate the normal to 
the surface from a point. 

Provides a smoothing algorithm for a tabulated function. 

Initializes an array to a specified value. 

Performs a binary or buffered read from the specified I/O unit into an 
array. 

Performs a binary or buffered write of an array onto the specified I/O 
unit. 

Rewinds the specified I/O unit. 

Computes the sum of all the numbers in an array. 

Reads a matrix stored by rows on the specified I/O unit into core, storing 
it in the standard FORTRAN order. 

Reads a specified number of rows of a matrix stored on the specified I/O 
unit into core, storing it by rows. 
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SYMBOLS 


Aero 

Symbol ; 

Data 

Group 

Data 

Item 

a 


ANGL 


ANGL 

NANGL 



• lANGL : 

ct 


CTOT 

c 

CHRD 

CHRD 

NCRD 

Cp ! 


CPRS s 


CPRS 

— — ■' 



NCPR 

Cf 


CSKN 


CSKN 

NSKN 

5 

DELB 

delF 



NDLB 




8* 

DELS 

DELS 

NDLS 



Angle of attack in degrees 
Number of angles of attack to 



Index of ANGL of present angle of 
attack 


' Total airfoil chord length 
Component chord length 
Index of point of maximum chord 

Surface pressure coefficient 
Index array for CPRS { 

Skin friction coefficient 
Index array for CSKN_ 

Boundary layer thickness 
Index array for DELB 

Boundary layer displacement 
thickness 

Index array for DELS 



























Aero 

Symbol 


Data 

Group 



FSMACH 


NMACH 



Freestream Mach number 

Number of freestream Mach , 
numbers to process 

Index of FSMACH of present 
freestream Mach number 


Boundary layer shape factor 


Index array for HSHP I 


Local Mach number 


Index array for ML 


Index of main component 

Index array of components being 
lofted 

Index array of reference component 

Index array of pivot points for lofted 
component 

Index array of pivot points for 
reference component 

X coordinates of pivot points 

Z coordinates of pivot points 

Flag to indicate lofted components 

Angle of rotation for lofted 
component relative to reference 
component 
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Aero 

Symbol 






Reft 10 

Pr 

k 

Po 


Data 

Data 

Group 

Item 


DRAG 


LIFT 

FARMS 

PITCH 1 


AXIAL 

j 


NORML 


SF i 

SCALE 

CREF 


SLOC 1 

SLOC 

.NSLC _ 


SLOT 

SLOT 

- — 


NSLT 


SRC’ 

SRC - 

NSRC 


TO 

SRFX 

RN 


PR 


KF 


Pp 


Explanation 


Drag coefficient 
Total lift coefficient 
Pitching moment coefficient 
Axial force coefficient 
Normal force coefficient 


Scale factor to convert input 
geometry data to feet 

Reference chord in feet 


Arc length (surface fitted 
coordinate) 

Index array for SLOG 


Location and height of slot exit for 
each component 


Index of first surface point beyond 
slot exit 


Source strength 
Index array for SRC 


Freestream stagnation temperature 

Reynolds number in million/feet 
Prandtl number 
Heat transfer factor 
Stagnation pressure 
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Data 

Group 

Data 

Item 

Explanation 


XSTAG 

X coordinate of stagnation point 


ZSTAG 

Z coordinate of stagnation point 

STAG 

ISTAG 

Index of stagnation points 


SSTAG 

Arc length of stagnation points 

THETA 

THETA 

Boundary layer momentum 
thickness 


NTHA 

Index array for THETA 

TRANS 

UGTRN 

User to global coordinate 
transformation data 


GLTRN 

Global to local coordinate 
transformation data 


XTRAN 

X coordinate of transition point 


ZTRAN 

Z coordinate of transition point 

TRNSIT 

LTRAN 

Fixed transition flag 


SSTRN 

Arc length of transition point 


VCMP 

Compressible surface velocity 

VCMP 

NVCM 

Index array for VCMP 


VSRF 

Incompressible surface velocity 
(corner points) 


NSRF 

Index array for VSRF 

VSRF 




VSRX 

Tangential component of 
incompressible surface velocity 
(corner points) 


NSRX 

Index array for VSRX 


























Aero 

Symbol 

Data 

Group 

Data 

Item 

i 

Explanation 

Vn 


VSRZ 

Normal component of incompressible 
surface velocity (corner points) 



NSRZ 

Index array for VSRZ 





Xg 

XGBL 

XGBL 

X coordinate, global coordinate 
system 



NXGB 

Index array for XGBL 

Xl 


XLOC 

X coordinate, local coordinate 

XLOC 

system 



NXLC 

Index array for XLOC 

, 

XUSR 

XUSR 

X coordinate, user input coordinate 
system 



NXUS^ 

^ Index array for XUSR 

Zg ' 

ZGBL 

ZGBL 

NZGB 

Z coordinate, global coordinate system 

r*" 

Index array for ZGBL 

'z.- 

ZLOC 

ZLOC 

NZLC 

i Z coordinate, local coordinate system 
Index array for ZLOC 

Z, 

ZUSR 

ZUSR 

Z coordinate, user input coordinate 
system 



NZUS i 

Index array for ZUSR 
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C/D > 


ero 

ymbol 


NCJX 

CJZ 


Explanation 


Source stream function AIC work 
array 

Index array for ASS 


Source X velocity AIC work array 

1 

Index array for ASX 

1 

Source Z velocity AIC work array 


Index array for ASS 


Vortex stream function AIC work 


array 


Index array for A VS 

i 

Vortex X velocity AIC work array 


Index array for AVX 

i 

Vortex Z velocity AIC work array 

j 

Index array for AVZ 

.1 

^ 




Work array for row of Jacobian 
matrix for wake centerline update 

Index array for AIJ 

Wake centerline chord lengths 

Index array for CJJ 

Wake centerline segment length in 
X direction 

Index array for CJX 

Wake centerline segment length in 
Z direction 


Index array for CJZ 
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Aero 

Symbol 


Explanation 


Convergence control flag 


Wake centerline segment angles 

Index array for CTH 

Stream function value array for 
wake centerline update 

Index array for STR 


Dynamic storage reference array 

Number of words in DYNM 
available for dynamic storage 

First available word in DYNM for 
dynamic storage 


Display area for error message 


LOCERR ’ Local error flag 

GLBERR ! Global error flag 


ETAL I I ; Trailing edge angle 


IXGB 

IZGB 


IVCM 

ILMH 

ICPR 

IDLB 

IDLS 

ISHP 

ITHA 


I Index for XGBL 
j Index for ZGBL 
1 Index for SLOC 

j 

I Index for VCMP 
Index for ML 
Index for CPRS 
Index for DELB 
Index for DELS 
Index for HSHP 
Index for THETA 
Index for CSKN 
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< c? 


ero 
ymbol 


Data 

Group 

Data 

Item 

Explanation 


lOTAB 

Table of available I/O units 

lOTAB 




NPRT 

Print unit 


ITNUM 

Number of current iteration 

ITNUM 




ITMAX 

Maximum number of iterations 


KEYA 

Process AVS, AVX, AVZ, ASS, 
ASX, ASZ t5T)e AIC’s on unit 1 


KEYB 

Process AVS, ASS type AIC’s in 



unit 1 


KEYC 

Process ASS, ASX, ASZ type AIC’s 

KEYS 

on unit 1 


KEYD 

Process ASS type AIC on unit 1 


KEYE 

Process AVX, AVZ, ASX, ASZ type 



AIC’s on unit 2 


KEYF 

Process ASX, ASZ type AIC’s on 
unit 3 

LEND 

LEND 

End of user input flag 

NC 

NC j 

Number of components 


NSP 

Maximum number of surface points 


NSM 

Minimum number of points on a 1 



surface of a component 

NSP 

NPN 

Number of surface points defined by 



user for each component 


IDISP 

Flag for input points redistribution 

NCMP 

NCMP 

Number of computational points 


28 






























Data 

Group 



Explanation 


Right-hand side array for vortex 
strength calculation 

Index array for RHS 


Work array for rows of the solution 
matrix 


Index array for SOL 


VARIN 


VARIN 


Work array for integration variable 
in wake and modified Goradia 
solution ‘ 


ICMP 


Number of component for which 
viscous flow solution is being 
performed 


ITRN 


ICORE 


Surface on the component 
Transition to turbulent flow flag 
End of potential flow core flag 


XWKA 


NXWA 


X coordinate, work array ] 
Index array for XWKA , 


XWKB 


NXWB 


X coordinate, work array 
Index array for XWKB 


XWKC 


X coordinate, work array 


NXWC 


Index array for XWKC 




















Aero 
Symbol ' 

Data 

Group 

Data 

Item 

Explanation 



ZWKA 

Z coordinate, work array 

Za 

ZWKA 

NZWA 

Index array for ZWKA 



ZWKB 

Z coordinate, work array 

Zb 

ZWKB 

NZWB 

Index array for ZWKB 

Zc 

ZWKC 

ZWKC 

Z coordinate, work array 



NZWC 

Index array for ZWKC 



DWAK 

Width of wake at the end of the core 



UWAK 

Wake velocity at the end of the core 

Xe ! 

WAKE A 

XWAK 

X coordinate of the end of core 



NWAK 

Index of surface point on aft 
i component at end of core 



IWAK ^ 

Number of wake centerline points 
before end of core 



















PROCESSING OF USER INPUT (A) 


In this section, the input cards and the preprocessing of the geometry data are 
described. The computer code is contained in OVERLAY (1,0), subroutines 


INPTR 

SMOPT 

READIT 

GEOMB 

GEOM 

LOFTR 

GEOMA 

SLOTR 

RESCL 

GLOBD 

AFPRM 

SRFIT 

COMPT 

GEOMC 


ANGLE WAKJlCt- 

LOCLD 

SURFD 

TRANE 

DISTP 

FTLUD 

DIR 


INPUT CARDS (AA) 


The input cards are read by subroutines INPTR and READIT. 


V Card 1 
* Title 

Card 2 

NC 

NSP 


DESCRIPTION OF INPUT CARDS 

Format (8A10) 

- 80 column title 

Format (2F10.0) 

- Number of components (1<NC<10) 

- Number of computational surface points, 21NC ^^ NSP ^ 200 


Cards 3 through 6 are input NC times 


Card 3 Format (2F10.0) 

NPP — Number of pivot points connected to this component 

(up to 3 pivot points per component) 

NPT - Number of points to be input for this component (sum of NPT 

for all components^306) 

Card 4 

(Xp,Zp) 


Format (6F10.0) 

- Coordinates of the pivot point referenced to this coordinate 
system; if NPP=0, skip this card 
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Card 5 
FMT 


Card 6 

(Xi.Zi) 


Card 7 
IM 


Format (8A10) 

- The format of the input point coordinates, enclosed in 
parentheses. Example (6F10.0) 

Format from Card 5 

- The input point coordinates in (Xi,Zi) pairs starting at the 
upper-surface trailing edge and ending at the lower-surface 
trailing edge 

Use as many cards as necessary 
Format (FIO.O) 

- Index of the main component 
Card 7 is skipped if NC-1 


Card 8 


Format (5F10.0) 


IC 

IPP 

ICR 

IPPR 


DELTA 


- Index of this component 

- Index of the pivot point to be used in placing this component 

- Index of the reference component 

- Index of the pivot point on the reference component to be 
used in placing this component 

“ Angle of rotation between this coordinate system and the 
reference coordinate system in degrees 


Note: 


Card 9 
NA 


This card is included for each component other than the main 
component, i.e., Card 8 is repeated (NC-1) times. Card 8 is 
skipped if NC=1 

Format (FIO.O) 

- Number of angles of attack to be input. (1=^NA^10) 


Card 10 

ALPHA 


Format (5F10.0) 

- Angle of attack in degrees, NA values 


Card 1 1 


Format (FIO.O) 


NM 


- Number of freestream Mach numbers. (I^NM^IO) 


Card 12 

FSMCH 


Format (5F10.0) 

- Freestream Mach number, NM values 


Card 13 


Format (2F10.0) 


CREF - Reference chord in feet. This number is used to 

nondimensionalize output quantities. 

SF - Scale factor. This factor converts the input geometry to feet. 
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Card 14 

Format (4F10.0) 


TO 

o 

- Stagnation temperature - R 


RN 

- Reynolds number - millions/ft 


PR 

- Prandtl number (use 0.77) 


KF 

~ Heat transfer factor (use 1.0) 


Card ir> 

Format (3F10.0) 


LTRAN 

- Fixed transition option, 

^ 0. r 1 implies free transition 
= 1. 1 jimplies fixed transition 


^XTRAN. 

' ZTRAX 1 

- location of fixed transition (use (0.0) if free transition' 


Card 15 is repeated (2NCi times. Upper 

suriace. firs: 


component: lower surfacte first component: 
second component: etc. 

upptT surfac(\ 

(’ard 16 

Format lAlOi 


THEbENDbbb 

“ The last data card of last case to be processed 
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SYMBOLS OF INPUT CARDS 


The following list of symbols is included to facilitate cross references to other sections of 
the computer code. 


r Theory 
► 

Code 

Definition 

^ef 

CREF 

Reference chord in feet 


IC 

Indices of components in the order that their data are 
stored 


ICR 

Index of reference component for each component \ 


IM 

Index of main component | 


IPP 

Index of pivot point used in placing each component 


IPPR 

Index of pivot point on reference component to be 
used in placing each component 

k 

KF 

Heat transfer factor 


^ LTRl^N 

Transition option: = 'O free Transition, = 1 fixed j 

transition 1 

1 Moc 

FSMACH 

Freestream Mach number 

1 

NA 

Number of angles of attack 

Nc 

NC 

Number of components 


NM 

Number of Mach numbers 


NPP 

Number of pivot points for each component 

• 

r 

% 

NPT 

Number of input points for each component 

I 

NSP 

Total number of computational surface points 

( Pr 

PR 

Prandtl number i 

1 Reft lO’** 

RN 

Reynolds number in millions per foot [ 


SF 

Scale factor for conversion of input geometry to feet j 
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Theory 

Code 

Definition 

: T„ 

TO 

Freestream stagnation temperature in °R 

j 

X,, Z, 

X, Z 

Surface point in input coordinate system 

Xp, Zp 

XP, ZP 

Pivot point coordinates in input coordinate system 

Xtr, Ztr 

XTRAN, 

Location of fixed transition in input coordinates 


ZTRAN 


O' 

ALPHA 

Angle of attack in degrees 

A 

t 

DELTA 

Angle of rotation from the component’s coordinate 
system into its reference component coordinate 
system ^ ,,, 
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PREPROCESSING OF GEOMETRY DATA (AB) 


The geometry package of the program is described in this section. The code is contained 
in the following subroutines. 


GEOM 

SMOPT 

SURFD 

TRANF 

GEOMA 

GEOMB 

GEOMC 

FTLUD 

RESCL 

LOFTR 

ANGLR 

DIR 

AFPRM 

SLOTR 

LOCLD 

SRFIT 

COMPT 

GLOBD 

DISTP 

WAKECL 


GEOMETRY 

DEFINITION 



j uf mult leh^ment airfoil geometries are contained in figure 7. Each of the 

j aiiioil gfonielnes shown rt pi'csents a two-dimensional cut thrtmgh a high lift wing 
configurat ion and ma\ consist of up to 10 airfoil components. The aa ioil geometries 
may be c|uite general having 

• Arbitrary distributions of camber and thickness 
0 Hiunt or pointed trailing edge shape 
I • Positi\ e. YA'Yo. or nc‘galiv( overlap of neighboring airioil aaou.'- 


f igure s illustrates some of these geometric features. ■ 

Note: Even though the geometry package of the program can handle 

quite arbitrary geometries, the user of the program must be 
warned that severe limitations are imposed on the geometry by 
the various aerodynamic models of the method. As an example, ^ 

the assumption of attached flow requires a smooth geometry ' | 

without abrupt changes of the airfoiTs surface. Other i 

limitations are pointed out in the description of the | 

ae rodynamic theory. 

The input geometry of an airfoil is defined by a set of surface points (Xj, Zj). The 
coordinates of these points may be specified in a different coordinate system for each 
j component. They are read into the program beginning at the upper surface trailing edge 

! point. The reading of the data then proceeds along the upper and lower surfaces of the j 

j airfoil component and ends at the lower surface trailing edge point. The trailing edge 
point of airfoils with a sharp trailing edge, appears twice in the data set. ' 

COORDINATE SYSTEMS 

The program uses three types of coordinate systems (fig.9). They are defined as follows. 

• Input coordinates (Xi,Zi) 

{ These are Cartesian coordinate systems selected by the user. The user can either ‘ 
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Figure 7. — Examples of Multielement Airfoil Geometries 



(a) Trailing Edge Closure 



Blunt 


(b) Overlap . 



Positive overlap Negative overlap 


Figure 8. — Geometry Features 


specify a different coordinate system for each component or can choose to define all 
input geometries in one and the same coordinate system. 


• Global coordinates (X<^,Zq) 

This is the input coordinate system of the main component. 

• Local coordinates 

There are two types of coordinate systems (fig. 9). One type is the boundary layer 
coordinate system, and the other is the local Cartesian coordinate system. 

The coordinates (Xj,Zj ) are defined such that their origin is located at the leading 
edge and the Xj axis passes through the average trailing edge point. The leading 
edge point is that point of the airfoil which is farthest away from the average 
trailing edge point. The distance between the two points defines the chord length, 
c, of the airfoil component. 

The coordinates (x,y) are surface fitted or boundary layer coordinates. The 
X coordinate is identical with the arc length, s, calculated by summation of the 
distances between surface points. The calculation of arc length begins at the lower 
surface trailing edge point. 
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Computational Surface Points 


To obtain accurate potential flow results, computational surface points are chosen which 
differ from the input surface points in both number and location. The calculations are 
performed in subroutine COMPT. 

The total number of computational surface points NSP is an input variable. The 
numbers of computational surface points N[ of an airfoil component is calculated using 
the formula 

N-(nSP- 2I Nc)|j; + 2I ,, 

where c is the chord length of the considered airfoil component and c^ is the sum of the 
chord lengths of all components. The symbol denotes the total number of airfoil 
components. The numerical result of the term 

(NSP- 21 

is truncated to its integer value. The above formula divides the NSP computational 
surface points among the Nq airfoil components such that each component has an odd 
number of points with a minimum of 21. 

The location of the computational surface points is determined in subroutine DISTP 
based on surface curvature. The new surface .points cluster in the regions of high 
curvature. The following sequence of calculations is used for each airfoil surface. 


Step 1 

The input surface points (Xj, Zj) j = 1, 2,... jmax used to define the arc length from 
the trailing edge , . ^ 

Si = 0 

^ ' (2i 

p 2 , ,2-1 V2 

Step 2 

The arc length is considered as a parameterization of the (X,Z)-coordinates of the 
surface X + X(s); Z Z(s). The curvature of the point (Xj, Zj) = (X(Sj), Z(sj)) is 
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where the derivatives on the right hand side are evaluated at Sj. 


Step 3 


A function SUM(K(s)) is computed using the formulti 


S '/4 

SUM(s)= f [K(X)] dX 

0 


(4 


Step 4 


The arc length is considered as a iuriction of ST M. Evaluation points are chosen as 
follows. ^ SUM (s- ): 




SUM: j=l,2,...J 

L.^ 


max 


(5 


Interpolation is used to predict :he arc length at SUM'^'j, 

Sj* = s(sUMj*) 


Note: 


If SUM represents the integral of curvature along the surface, 
then \ are separated by equal increments of curvature. 


Step 5 


Finally, interpolation is used to compute the values of X,Z corresponding to s'" j. i.e.. 


Lofting 



Jmax 


The components other than the main component are rotated and translated from their 
input coordinate systems into the global coordinate system. This is done in subroutine 
LOFTR as follows. 


For each component a reference component is specified. Translation is performed by 
moving a specified pivot point in the coordinate system of the translated component to a 
specified pivot point in the coordinate system of the reference component. Rotation is 
performed about the pivot point by a specified angle between the reference component 
coordinate system and the coordinate system of the rotated component. 

The transformation is performed according to these equations: 


Xj = ^Xj ~ Xp j cos A + (Zi - Zp) sin A + Xp 
Z j — ^ Zj — Zp^ cos A — ^ Xj — Xp^ sin A + Zp 



The symbols have the meaning: 

Xj, Zj Surface point in coordinates of reference component 
Xj, Zj Surface point in coordinates of lofted component 
A Angle of rotation ' posit ivo clockwise ^ 

Xp,Zp Pivot point in coordinates of lofted component 
X^, ZJ, Pi vot point in c(U)rdinates of reference component 


The program ensures that a component geometry is not rotated and traiislatec. in: :h- 

global coordinate system unless its reference component has been rotated and 
translated. 

As an example, figure 10 shows the lofting of a high lift airfoil with four compoi ents. 
The second component is the main component (wing). Its coordinate system serves as 
the global coordinate system in which the geometry of the other three compi>r;ents 
(leading edge flap and trailing edge flaps) has to be defined by the process of ltdfing. 
The points A, B, C are the pivot points. The rotation angles are also indicated in the 
figure. The following table lists the components and their reference components The 
sequence of the lofting procedure is the sequence in which the components arc shown 
below. 


Component 

Reference Component 

A 

Pivot Point 

1 

2 

1-2 


4 

L 

4-2 

B 

3 

4 

3-4 

C 


The pivot point, C, is transformed into the global coordinate system during th - If ling 
of component 4. 


AIRFOIL PARAMETERS 
Slot Height 

The slot height at the exit of the slot is defined as illustrated in figure I 1 lor the two 
geometric cases of positive and negative overlap of neighboring airfoil component>. A 
straight line defining the slot height is drawn from the lower surface trailing edge point 
of the upstream component perpendicular to the surface of the downstream component. 

The slot height, hsiot' ^^id the coordinates of the point, P\, on the surface of the aft 
component are determined in subroutine SLOTR as follows: 
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+ A 

Component 1 





Component 2 
(main component) 





Component 4 
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Step 1 


The minimum distance Ptf?! from the trailing edge point, Pry:, to the points on the 
upper surface of the aft component determines the point Pj (fig. 12). 



Figure 12. — On the Calculation of Slot Height 

Step 2 

A curve, Cj, passing through the points Pi— i, Pi, Pi+i is calculated. The parametric 
representation of Q in terms of a nondimensional chord length ^ reads 

Xi(r) = rXi + (i-n Xj^i +f(i -d aj 

( 7 ) 

Zj (f) = rZi + (i -?) Zi_i + f(i 

The values of ^ at the various points are shown in figure 12. The point, Pg., w^here 
^ = aj,is the projection of Pi+i on the chord length The variable aj is obtained from 
geometric relations as 
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Since the computational surface points are distinct, one always computes aj>^ 1. The 
coefficients aj, /3j of the equations for Cj follow from 


_ Xi-n-aiXj-(l-ai) X-_, 
“i aj (1 - aj) 


/3i = 


_^i+i 0 


aj (1 - aj) 


Step 3 


The coordinates of the point P\ = (X\\ Z\) are calculated using 

Zn " fN^i + - ^n) ^i-1 + (* - 


(9) 


(10) 


and the equation of the normal to the surface at that passes through the 
trailing edge point (X^iv K t 



Combining these equations produces a cubic equation for 

Cq+ C j + 

I 

— - ■ ■' .. - III II ■fP " WPIIfP H I ■ ^ IIP ■ I ■ 

with the coefficients 

c, = - 2»i(XTE - Xi_, ) - 2(Ji(ZxE - Zj., ) - ( X| - Xj_, + «i)2 - ( Z, - Z,., + J,)2 

C2 = 3aj(Xj-Xj_j +0!i)+3^i(Zj-Zi_i +^j) 

C3 - - 2 (aj- + ) 

The real value of in the range 0 ^ =saj is chosen. Once Cn is known, the 
coordinates can be calculated from equation (10). 

Step 4 

The slot height follows from 

^slot" W^TE"%) +(^TE"^n) (12) 
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Trailing Edge Closure Angle 


Given the coordinates of the computational surface points Pj, P2, Pns Pn+i iri global 
coordinates, the trailing edge closure angle of the m>th airfoil component is 
calculated from ^ 5- 

'm = "C '131 


where Sj, Sj^are the slopes of the first and last segments of the airfoil component, i.e., 


c 

dZ 

i 

^1 

dX 

X 

1 

1 


dZ 

„ Zn+1 - Z]SJ 

t>N 

dX 

N ^N +1 " % 


The notation is illustrated in figure 13 . 


Zg 

k 



Xr 


Figure 13. — Airfoil Trailing Edge Geometry 



AERODYNAMIC ANALYSIS 


The flow field of multielement airfoils can be divided into several flow regions with 
j different physical characteristics. These are: 

■ m* Outer potential flow 

• Laminar boundary layers 
r 

■ '• Transition 

* % Ordinary turbulent boundary layers 

I 

\ • Wakes 

k 

k 

[ • Confluent boundary layers 

# Regions of separated flow. 

The last two flow regions distinguish the flow problem of multielement airfoils from the 
problem of airfoils at cruise conditions. 

The following sections of this document contain a detailed description of the 
mathematical formulation and solution of each flow region. Regions of separated flow 
are not modeled by the aerodynamic analysis. 



ITERATION PROCEDURE (BA, BB, BC) 


GENERAL DESCRIPTION 

The concept of displacement thickness is used to represent the effect of the various 
viscous layers on the outer potential flow. Instead of adding the displacement thickness 
to the airfoil geometry, a distribution of sources along the airfoil surface and along the 
wake centerlines is utilized for the simulation of the viscous flow displacement effects. 
This is the so-called surface transpiration method which, within the framework of thin 
boundary layer theory, is completely equivalent to the method of adding geometrically 
the displacement thickness to the basic airfoil geometry. Details of the scheme are 
described in the section titled Viscous Flow Representation. 

The mathematical formulations of the inviscid and viscous flow problems are coupled 
through their respective boundary conditions. A solution of the potential flow problem, 
such as the distribution of surface velocity, depends on the airfoil geometry and on the 
boundary layer displacement thickness. The solutions of the boundary layer problems, 
in turn, which include the displacement thickness, depend on the potential flow 
velocities. 

The objective of the solution procedure, therefore, is to find those particular 
distributions of surface velocity and boundary layer displacement thickness which 
simultaneously satisfy both the potential flow problem and the boundary layer 
problems. The desired solutions of surface velocity and displacement thickness, from 
which all other flow parameters can be computed, must be arrived at in an iterative 
procedure since the coupling of the flow problems is mathematically nonlinear. The 
computer program uses a cyclic iteration procedure described below. The main loop of 
the iteration procedure is contained in subroutine SYSOL. 

CYCLED 

The computation performed during this iteration cycle consists of the following steps. 

1 . The first potential flow solution is calculated without any representation of viscous 
flow effects. 

2. The position of the wake centerlines is computed. 

3. Solutions of all viscous flow problems including laminar and turbulent boundary 
layers, confluent boundary layers, and viscous wakes are calculated with the 
surface velocities and wake centerline location obtained in the previous step as 
input data. At the end of this computational step, a first estimate of the 
displacement thicknesses' of all boundary layers and wakes is available. 
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CYCLE 1 


The following computational steps are performed. 

1. A source distribution representing the displacement effect of all boundary layers 
and wakes is calculated. 

2. A new potential flow solution is calculated for the basic airfoil geometry with a 
distribution of sources computed in the previous step along the surface and wake 
centerlines of the multielement airfoil. 

3. The position of the wake centerlines is updated using the result of the previous 
potential flow calculation. 

4. All boundary layer and wake properties are recomputed with the last available 
potential flow velocities and wake centerline geometries as input data. 

In subsequent cycles of the iteration procedure, the calculations described under Cycle 1 

are repeated. Several refinements of the iteration procedure which are used to assist in 

the convergence of the scheme are described in the following chapter. 


CONVERGENCE ASSISTANCE 

The following techniques are used to improve the convergence characteristics of the 
iteration procedure. 

SMOOTHING 

The distributions of surface velocity and displacement thickness exhibit rapid changes 
or discontinuities in certain areas, due to deficiencies of the chosen aerodynamic model, 
and are physically not realistic. Anomalies of this kind typically occur in the following 
areas. 

• Trailing edges of airfoil components 

• Transition points 

• End of the potential core region 

Consequently, the source distribution, which is computed from the potential flow 
velocity and the boundary layer displacement thickness, will also exhibit an unrealistic 
behavior in these areas. To arrive at a physically meaningful converged solution, the 
computed source distributions are modified as described in detail in the section titled 
Viscous Flow Representation. 



SCALING 


To assist the iteration scheme to arrive at a converged solution, the increments of the 
computed source strength cr are scaled according to the formula 

(j(i) = 

The symbol denotes the computed source strength of the i-th iteration cycle prior to 
scaling. The scaled source strengths of the i-th and (i-l )-st iteration cycle are denoted by 
and respectively. 

The formula states that the source strength a is scaled by adding 2/3 of a computed in 
the present iteration cycle to 1/3 of a a computed in the previous iteration cycle. 

CONVERGENCE CHECK 

The code does not rely on a convergence criterion. Instead, all solutions are obtained in 
five iteration cycles whether or not a converged solution is arrived at after the last 
cycle. The quality of the convergence of the solution must be judged by the user of the 
code. 


A reliable check on the convergence of the iteration procedure does not seem to exist. A 
convergence check on lift coefficient and/or surface pressures might be useful for some 
cases but is not always reliable. An automatic check on lift coeffient could terminate 
the computation before a truly converged solution is obtained. Nevertheless, the 
computer code contains a dummy subroutine called CONVR which can be used for the 
addition of a convergence check. 
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VISCOUS FLOW REPRESENTATION (BBA) 


' The representation of the displacement effect of boundary layers and wakes in the 
solution procedure for the inviscid part of the flow field is described in this section. The 
computer code is contained in OVERLAY (0,0), subroutine SRCER. 


EQUIVALENT SOURCES 


The method is the so-called surface transpiration method in which a distribution of 
sources along the airfoil surface and along the wake centerlines simulates the 
displacement effect of the viscous layers. The strength o- of this equivalent source 
distribution is calculated from 



(15) 


In this formula, cr denotes an incompressible source, nondimensionalized by the 
freestream velocity Uoo- The symbol Vj stands for the incompressible dimensional value 
' of the surface velocity. 


The use of the boundary layer displacement thickness in the equation for the 
computation of a requires a detailed explanation. In the computer code, the 
displacement thickness 8* is calculated using a mixed compressible-incompressible 
method. This approach is adequate within the theoretical framework of the 
Karman-Tsien compressibility correction which does not require any geometry scaling 
and, consequently, does not distinguish between compressible and incompressible values 
of the displacement thicknesses. 


MODIFICATIONS OF SOURCES 

The computed distribution of 8* is discontinuous at certain locations, e.g,, at the point of 
transition from laminar to turbulent boundary layer flow, and at the end of the core region. 
Furthermore, the distributions of displacement thickness and potential flow velocity exhibit 
a rapid change near the trailing edge of an airfoil. These discontinuities and rapid changes of 
8* and Vj are caused by deficiencies of the aerodynamic modeling and are not physically 
realistic. Consequently, the source distribution, which is computed from 8* and Vj will also 
exhibit such anomalies in certain areas, and would produce erroneous results or could even 
lead to catastrophic program failures if left uncorrected. Therefore, the following 
modifications of the source distribution are made. 


Before the sources are computed, the distribution of 8* on all airfoil surfaces is 
smoothed once. Wake displacement thicknesses are not smoothed. 


A negative value of an airfoil source is eliminated by substituting the last positive 
source strength found upstream, i.e., in the direction of the stagnation point. 
Negative wake sources (sinks) are not modified. 
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• Unlimited growth of the source strength, upstream and downstream of a trailing 
edge, is avoided by overriding the computed values of cr. Airfoil sources are kept 
constant on the last four segments of each airfoil surface. Similarly, wake sources 
are constant on the first four segments of the wake centerline. The distributions of 
airfoil sources and wake sources are continuous, but cr will in general be 
discontinuous at trailing edges. 

• The value of the source strength cr is limited to the range - .07 ^ cr ^ .07 
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POTENTIAL FLOW (BBB) 


The calculation of the incompressible potential flow solution, the compressibility 
corrections, and the calculation of the stagnation points are described in this section of 
the document. The computer code is contained in OVERLAY (3,0) consisting of the 
following subroutines: 


POTLF 

POTLFE 

AICVAL 

WAKET 

REDUCE 

POTLFA 

ANLYS 

CPTAIC 

WAKEJ 

WAKES 

BCKSUB 

QNEWT 

POTLFB 

COMPR 

SAVAIC 

WAKEG 

NEWTR 

POTLFC 

STAGN 

GETAIC 

NEWTR 

VIP 

POTLFD 

— - - 

GAUSSR 

RSEITE 

RECVEC 

EVAL 


METHOD OF OELLER 


The potential flow solution uses the stream function method of Oeller (ref. 11). Its main 
assumptions are that the flow is 

• Two-dimensional 

• Incompressible 

• Irrotational 

• Attached to the airfoil surface 


The principles of the potential flow method are now introduced using a single airfoil 
without boundary layer representation as an example. The problem is formulated and 
solved in global coordinates, (Xg,Zq), where the subscript G is dropped for convenience. 
The above assumptions allow the problem to be formulated in terms of the stream 
function as the dependent variable. The stream function ''T is governed by Laplace s 
equation. 


v-^'p= 0 


;(i6) 


which is linear, so the solution of the flow field can be obtained by superposition. The 
airfoil is represented by a distribution of vorticity along the surface of ^strength y(s). 
Adding the stream function of a uniform freestream, whose velocity Uoo meets the 
X(; axis under an angle of attack a, to the stream function of this vortex sheet results in 
the stream function of the whole flow field. 

1 ®TE 

'!' = UooCOsaZ-UooSino;X + ^ f 7 (s')Cnr(s, s')ds' ' (17) 

o"' ^ 

-V / ,, 

free stream vor^x_sheet ^ j 



The notation is illustrated in figure 14, where, in particular, the radius r is the distance 
between a point on the airfoil surface and a field point (X,Z). The stream function of the 
vortex sheet is found by integration of the stream functions of elementary vortices from 
the lower surface trailing edge point (s=0) to the upper surface trailing edge point 
(s= ste)- The value of the stream function 'k is constant along a streamline. Hence, 4^ 
is also constant along the airfoil’s contour, which is part of the stagnation streamline. 
This fact is used in calculating the unknown strength of the vortex sheet, y, and the 
unknown value of 'k at the airfoil surface from equation (17) 


To solve this integral equation, the airfoil geometry and the vortex distribution are 
discretized as follows (fig. 15). The airfoil surface is divided into N segments. The (N + 1) 
corner points of these segments are placed on the airfoil surface and are then connected 
by straight lines, i.e., the airfoil geometry is represented by a polygon. The vorticity is 
distributed along this polygon such that its value is constant along each segment. If 
further N collocation points or control points (Xj,Zi) are chosen, the integral equation 
(17) reduces to a set of linear algebraic equations. 

I 

N V . ^ 

^^-2 K:; 7 j = UooCos (x Zj - Uoo sin a Xj (i-l,2,...N) (18) 

j=l 



Figure 14. — Notation of Stream Function Method 


Segment corner point 



Segment j with 
vortex sheet of 
constant strength 7j 


Control point i at center 
of segment 


Figure 15. — Discretization Used in Potential Flow Method 


Here, is the value of the stream function at the contour of the airfoil; and the 

Ky’s 

are the aerodynamic influence, coefficients of a constant strength vortex sheet, defined 

by 


/^^*<^"r(Si,s;)ds' (19) 

The set of linear equations (18) contains N+1 unknowns, i.e., N unknown vortex 
strengths yj plus one unknown value of the stream function 'kp, but provides only N 
equations. The missing equation is supplied hy the Kutta condition, which is formulated 
as 


7, +7^=0 


( 20 ) 


i This equation is not the classical Kutta condition, which, in potential flow past airfoils 
with a nonzero trailing edge angle, postulates the existence of a stagnation point at the 
trailing edge. Instead, equation (20) states the velocities at the upper and lower surface 
I of the trailing edge are equal, but not necessarily zero as would be the case at a 
j stagnation point. 

VORTEX STRENGTH AND STREAM FUNCTION 


The principles of Oeller’s potential flow method are explained in the previous section. 
Modifications of the method to include source distributions along the airfoil surface and 
the wake centerline for the simulation of boundary layer and wake displacement effects 
are described next. — i- » - ~ 
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SINGLE AIRFOIL 


Sources are distributed along the airfoil surface and along the wake centerline to 
represent the displacement effect of boundary layers and wake. The curvature of the 
wake is assumed to be negligible, so that vortices are only placed along the airfoil 
surface and not along the wake centerline. The stream function of the flow field must 
include the contribution of this source sheet. 

Accordingly, equation (17) is modified to 

1 ^TE 1 ^E 

^ = Uqo cos q: Z ~ Uqo sin ol X S T(s') r (s,s') ds^ + J a(s') </>(s,s^) ds' (21) 

o o 

Here, the symbol o-Cs) denotes the strength of the source distribution. The angle (p and 
the arc length sjr are illustrated in figure 16. 

To solve for the unknown strength of the vortex sheet, the potential flow problem is 
discretized as follows. The airfoil surface is subdivided into N segments. In addition, 
there are N\\’ segments representing the wake centerline. Each of the N segments, 
which approximate the airfoil surface, carries a constant strength source sheet crj and a 
constant distribution of vortex strength yy Only source sheets of constant strength are 
placed on the N\v wake segments. 

Furthermore, to solve equation (21), a streamline of known position must be chosen 
along which the stream function will have a constant value Knowing the position of 
this streamline, control points (Xi,Zj) on the streamline can then be chosen. In the 
absence of sources, the airfoil surface itself represents without any doubt such a 
streamline. The question now arises: does the airfoil surface remain a streamline when 
sources are distributed along the surface? This is obviously not the case, since the 
sources model the displacement thickness of the boundary layer which, when added to 
the geometric airfoil surface, produce the so-called displacement body. That is to say, 
the presence of the sources has shifted the stagnation streamline from the airfoil 
surface to the surface of the displacement body. Figure 17 shows the qualitative pattern 
of the streamlines in the vicinity of the airfoil surface. 

Nevertheless, control points that are located on the airfoil surface and, therefore, are 
fixed during the solution procedure are chosen for reasons of simplicity and 
computational efficiency. This choice can be justified as follows. 

The objective of the potential flow calculation is to provide a solution of Laplace’s 
equation with a known distribution of sources in the flow field simulating viscous tlow 
displacement effects. The strengths of these sources represent boundary cond tions 
prescribing the velocity normal to the airfoil surface and the wake centerline. When 
discretizing equation (21), one can make use of the fact that the source (7i on the i-th 
segment already satisfies the specified normal velocity boundary condition. Hence, the 
superposition of the effects of all remaining vortices and sources of the flow fit Id must 
satisfy the boundary condition of zero normal velocity at the i-th control point Writing 
this condition as 6^/8x = 0, where x is the local Cartesian coordinate tangent to the 



surface, one arrives at the boundary condition of a constant value of the stream function 
along the airfoil surface. This constant value of ^ is denoted by 

The reader is reminded the potential flow problem in the presence of sources is solved 
as if the stream function had a constant value along the surface. The superposition of 
all vortices and sources of the flow field including a[ will produce a stream function 
whose value changes along the airfoil surface. 


Discretizing equation (21) as described results in 



Figure 16. — Additional Notation for Source Distributions 
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N Y N+Nw _ 

- 2 Kjj 7 j = Uoo cos a Zj - Uoo sin a Xj + 2 Ky oj 

j=l j=l 


(i=l,2, ... N) (22) 


where the stream function influence coefficients of sources are defined by 



1 /j+1 , -V . , 


(23) 


According to the argument preceding equation (22) , this equation is a superposition of 
the effects of all vortices and sources of the flow field on the control point of the i-th 
segment, but does not include crj, whose influence is eliminated by formally setting 
K\\= 0. The control points (Xi,Z[) are the midpoints of the airfoil segments. 


KUTTA CONDITION 


When adding sources to the flow field the formulation of the Kutta condition must be 
modified. Requiring that the velocities on the upper- and lower-surface trailing edge are 
equal results in 


7l = (^N~^l) sine 


(24) 


The variable e denotes the trailing edge closure angle (fig. 18). The symbols 
represent the source strength of the first and last airfoil segment, respectively. The 
reader should note that the source strength ctn+i of the neighboring wake centerline 
does not enter this formulation of the Kutta condition. 


MULTIELEMENT AIRFOIL 


The formulation of the potential flow method of a single airfoil is readily extended to 
airfoils consisting of several components. Noting that a multielement airfoil with N^. 
components has stagnation streamlines with different values of the stream function, 
one arrives at 


V 


Nc Nn, 

2 2 Kjj 7j = cosa(Zi) -sina(Xi) 

m=l j=l " 


m 


Nr 

! ^ 

, m=l 


(N+Nw)^ 
2 

j=I 


KjjOj 


(25) 


(m=l, ... N^) 
0=1,2, ...N^) 
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Airfoil 



Wake 

centerline 


Figure 18. — Airfoil Trailing Edge 


The subscript m indicates the m-th airfoil component. In particular, the coordinates ot 
the Nm control points on the surface of the m-th component are denoted by (Xj lm* (Zj 
and the stream function value of the m-th component is Furthermore, in equation 
(25) the singularity strengths yj, crj, as well as the stream function value represent 
nondimensional quantities, referred to the magnitude of the freestream velocity Ux- 


Equation (25) represents a total of 


N-y — 


N( 

2 N 


m=l 


m 


equations, but There are Nj + Nq unknown vortices yj and stream functions 
additional equations are provided by the Kutta condition of each airfoil component, 
which reads 




(26) 


where is the trailing edge closure angle of the m-th airfoil component. 

INCOMPRESSIBLE SURFACE VELOCITY (BBBAl) 


The incompressible surface velocity is computed in subroutine POTLFE. Equations for 
the velocity components U, W in global coordinates are derived from 


Nc 

+ 2 2 K; 

m= 1 j= 1 


(N+Nw)^ 


ij 


7j + 


TO C 

2 2 K;; a, 

m= 1 j= 1 




^ (X,Z) = cos a Z - sin a X 


(27) 



where yj, o-j, are nondimensional quantities, referred toj^l^. 

From the equations defining the stream function 

U: 


one obtains 


=(^\ 
u<« \az/i 

Uoo lax A 


u 


= COS a + 

°° m=l j= 


NC N, , S 

2 2 ^ Tj + 2 2 — y 

az m=i j=i az 


( 28 ) 


(29) 


Nr- N V Np (N+N^) S 
C 0K - ^ "^'m 3K-- 

1 . V V -V- y y 

-rr=sino:- 2 2 —J T| - ^ 


W; 

U 


m=l j=l 


m=l j=l 


ax J 


The velocity influence coefficients 


V V S S 

aKjj aKjj aKjj aKjj 

ax ’3X ’iz 


are discussed in the next section. 

The incompressible surface velocity follows from 


Vj //Uj Y /Wj 

" no;;/ ^VOT/ 


(30) 


(31) 


In computing the velocity components Uj, Wj at control points on the airfoil surface, the 
source induced velocity normal to the surface is not included 


0. 


ax 


= 0 . 


az 


Therefore, the surface velocity Vi isitangential to the airfoil surface. 

AERODYNAMIC INFLUENCE COEFFICIENTS 


The aerodynamic influence coefficients are computed in subroutine CPTAIC. 



AM FUNCTION INFLUENCE COEFFICIENTS 

These influence coefficients are defined by equations (19) and (23) for a vortex segment 
of constant strength and a source segment of constant strength, respectively. 
Introducing local segment coordinates, see figure 19, the influence coefficients can be 
rewritten as 

o 


K,f = / arc tan ^ d£j 

O <1 



Figure 19. - Local Segment Coordinates ^.n 



Hence, 


47TKy =^iCn(^^+T,f)-(^j-Cj)en [fe-Cj)^ +7,^] 
- 2 Cj + 2i?j ^arc tan arc tan ^ ^ ^ 

47t Kjj = 7?i Cn + T?f ) - Vi Cn - cj^ + 77fj + tt cj 

/ \ 

-2 arc tan — + 2 (f: - cj j arc tan — 

‘ \ ‘ V 


' ( 34 ) 


( 35 ) 


The control point coordinates (^i,i 7 i) are expressed in terms of global coordinates as 
follows. 




J 

-(Xi-Xj)(Zj„-Zj)+(Xj^,-Xi)(Zi-Zj) 

".= Cj 


with 



Here, (Xj,Zj), (Xj4-j,Zj+j) are the segment corner points in global coordinates and (Xj,Zj) 
is the control point in the same coordinate system. 

The stream function influence coefficients and Kj® are calculated using the above 
equations. The latter coefficient is set to zero for source segments on the airfoil surface, 
i.e., 
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r 


TT- ^ -f 




VELOCITY INFLUENCE COEFFICIENTS 

it 

The derivatives of Ky and Ky with respect to the global coordinates X,Z (subscript G 
dropped for convenience), are termed velocity influence coefficients. They are expressed 
I in terms of local segment coordinates f ,17 by 


V 

V 

V 

3Kij_ 

9Kij 


9X 


9X 9t7 9X 

V 

V 


9Ky 

9Ky 

9^ 9Ky 977 

9Z ■ 


9Z 9ti 9Z 

S 

S 

S 

9Kij 

9Kij 

9$ . 9t? 

9X " 

9^ 

9X 9tj 9X 

S 

S 

S 

9Kjj 

9Kjj 9Kjj 

9Z 

9$ 

9Z 9t? 9Z 


Here, the coefficients of the Jacobian of the transformation from local to global 
coordinates are 


9^.Xi+l-Xj 

9^ ^j+1 


9X Cj 

az Cj 

(38) 


■9t7_ 

" - — - - 

9t? _^j+l 

9X c j 


9Z ~ c j 


The velocity influence coefficients in local segment coordinates read 


V S 
9Kjj ^ 3K;, 


dr] 


i' {«n (jf + qf ) - tn |^( li - j j 


V s 

9 r/ ~ 3 ^ 


= arc tan— - arc tan ^ 


( 39 ) 


( 40 ) 


In calculating the influence of a segment singularity distribution on its own control 
point, the vortex velocity influence coefficients are computed at 17 = 0 + and the source 
velocity influence coefficients are set to zero. 
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( 41 ) 


^ 3Kjf ^ 

9X ^'9Z 


STAGNATION POINTS! 


A stagnation point is a point on the airfoil surface where the velocity Vj is zero. The 
location of the N(^ stagnation points is determined in subroutine STAGN by marching 
along the surface"^f each airfoil component from the lower to the upper surface trailing 
edge. A change in sign of the surface velocity Vj and linear interpolation between 
segment corner point coordinates determine the location of a stagnation point. 

COMPRESSIBILITY CORRECTION' 


Compressible potential flow solutions are calculated in subroutine COMPR. 


The well known Karman-Tsien rule (ref. 12) is used, which relates compressible and 
incompressible quantities of the same airfoil geometry. A scaling of the geometry is 
therefore not performed. 


The compressible surface velocity, V^, is obtained from the corresponding 
incompressible velocity by means of the equation 



where the parameter X depends on the freestream Mach number 


(42) 



(43) 


The way equation (42) is written indicates all velocities are nondimensionalized by the 
freestream velocity Further, the local Mach number, M^, and the surface pressure 
coefficient, Cp, are given by the isentropic flow relations. ~ 


M 


e 


M. 


V. 






(44) 


2 

cp- — j 

7^00 



0.001 


(45) 
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T=l+l^Mi l-(^) ;T=1,4 


The code uses Bernoulli’s equation 


■'-(if 


to compute the surface pressure coefficient for 0. 

At a stagnation point, compressible variables are calculated from 


v, = o 


Me = 0 


1 


Dimensional compressible surface velocities, V^^., in ft/sec are obtained from 


V =(^)u 


U«,= 49.02 Moo 


1 + .2M, 


F— 1 

|_sec J 


The symbol Tq stands for the freestream stagnation temperature in °R. 
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LAMINAR BOUNDARY LAYER (BBCA) 


This section describes the calculation of compressible laminar boundary layer 
characteristics. The method is coded in OVERLAY (4,1 h subroutine LAMNA. 

COMPRESSIBLE LAMINAR BOUNDARY LAYER 

COHEN AND RESHOTKO METHOD 

The method used in LAMNA is based on the compressible analysis of Cohen and 
Reshotko (ref. 13) as modified by Goradia (ref. 1). Some results of the Pohlhausen 
method (ref. 14) for incompressible boundary layers are also utilized. The following 
assumptions are made: 

• The laminar boundary layer flow is steady and two-dimensional 

• Curvature effects are negligible 

• Dynamic viscosity is a linear function of temperature 

• Surface temperature is uniform 

• Air is a perfect gas 

Consequently, the equations governing compressible, steady, two-dimensional laminar 
boundary layer flow read 


|j-(pu)+|^(p») = 0 


(50) 


pu 


9u 

9x 


3u ^ dp ^ 9 / 9u \ 
8y dx ^ 9y / 


(51) 


P c. 


/ 3T 3T\_ dp^3/.3T\. / 3u ^ 

(^3x ^ 3y ) ”dx 3y\ 3y / ^ \3y j 


(52) 
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where the symbols denote 


Cp Specific heat at constant pressure 

k Heat transfer factor 

p Static pressure 

T Static temperature 

u,v Velocity components in x,y directions 

x,y Surface fitted coordinates 

fx Dynamic viscosity 

P Densitx 

Two additional equations are given by the perfect gas law 


p =pRT 


( 53 ) 


and the assumption 


iL 

Mq 



( 54 ) 


The subscript o denotes freestream stagnation values. The coefficient X is determined by 
matching the viscosity with the Sutherland value at the airfoil surface. Sutherland’s 
viscosity formula reads 


ju 

■TTo^Tn^su 


+ *^Su/Tl\ 

^o/ 


3 

T 


( 55 ) 


Ksu-198.6°R 

Hence, 


- T + Ksu V To 


( 56 ) 


The subscript w denotes values at the airfoil surface. The listed equations are solved 
subject to the following known boundary conditions. At the airfoil surface 



u(x, 0) = v(x, 0) = 0 


T(x,0) = T^ 


and at the outer edge of the boundary layer (y = 6), denoted by the subscript e. they 


u(x, 6) = 


T(x,6) = Tg 


The code uses the integral method of Cohen and Reshotko for the solution of the 
described flow problem. An outline of this method is given below. 

The compressible boundary layer problem is converted into an equivalent 
incompressible problem by means of the Stewartson transformation, which reads 


^ _ ££ 

” Po 


9^ pv 

9x 
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X= / dx 

° Po 


Y=!i f 


O O 


(58) 


U = 


3^ 

3Y 


3X 


(59) 


The symbols 'I' and a represent the stream function and the speed of sound, 
respectively. Quantities of the equivalent incompressible problem are denoted by capital 
letters, i.e., U, V velocity components; X, Y surface fitted coordinates. 

The Prandtl number 


Pr = 



and the specific heat Cp are assumed to be constant during this transformation. 


All integral form of the momentum equation governing the momentum thickness and ; iu“ 
d is)jiacement thickness 9' t|. can be derived for the transformed problem. Defining 
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( 61 ) 


where 



is an enthalpy term involving the local stagnation enthalpy, and the freestream 
stai^nation enthalpy, h^, , and the symbol A stands for the boundary layer thickness in the 
translormed space, the momentum integral equation takes the form 


dX 


U. 


dUe/ *\ 

dX (2^tr + 5tr) 




(62) 


The symbol Uq is the freestream stagnation value of the kinematic viscosity. 


Equation (62) can further be expressed in terms of two parameters, the shear parameter 
S , defined by 



(63) 


and the correlation number n, defined by 


0tr dU. 


n= - 


dX 


Also, introducing the shape parameter 


^tr = e 


tr 


yields the following momentum integral equation 


(64) 


e dx/duA 


= 2 


,dX / 


[n(Htr+2)+«] 


(65) 


This equation is the basic equation of the solution method for the compressible laminar 
boundary layer calculation in subroutine LAMNA. It is solved for n assuming that the 
RHS is a function of n and Sw only, i.e., 


N(n,Sw)-2[„(H„+2)+t] 


( 66 ) 
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and that this functional relation provides sufficient accuracy by similar solutions 
investigated in reference 15 by Cohen and Reshotko. 

The parameter S\v is the value of the enthalpy at the airfoil surface. 


S 


w 



(67) 


Numerical values of the function N(n,S^are listed in table 2 and are plotted in figure 
20 between the stagnation point and the separation point for two surface temperatures. 
The latter figure shows that the function N(n,S\v) can be approximated by 


N = A + Bn 


(68) 


for fixed values Sw- Hence, equation (65) can be integrated resulting in 

X 


n = * A U, 


_B d Ue f U® ^dXi 
dX 


(69) 


Finally, equation (69) can be expressed in terms of physical compressible variables by 
means of the Stewartson transformation, (eqs. (57), (58), (59)). 


n = - A M 


-B d Mg 
e dx 


_4 X B-1 4 

Te f Me dx, 
o 


(70) 


Here, Mg = Ug/ag is the local Mach number at the outer edge of the boundary layer and 
Tg is the temperature at the same location. Both Tg and Mg are related by 



(71) 


Ec)uation (70' is evaluated numerically in the computer code for known dislribuiions of M,.. 

Having computed the correlation number, n, the other boundary layer parameters of 

inii'iest i6‘ . w, C|, Rg ) can be calculated. Details of the calculations follow. 

% 

COMPUTATIONAL DETAILS 


Step 1 

The parameter S\v is obtained from 


^-1 1 2 
1 + ^^rMoo 


(72) 
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Sw 

n 

i 


- 0.8 

0.1215 

0 



0.1304 

0.0312 



0.1298 

0.0436 



0.1260 

0.0681 



0.1212 

0.0827 



0.1017 

0.1214 



0.0355 

0.1935 



0 

0.220 



- 0.0837 

0.2678 



- 0.2008 

0.3179 



- 0.2522 

0.3366 


- 0.4 

0.0899 

0 




0.0300 




0.0624 




0.1210 



0 

0.220 



- 0.0722 

0.3019 



- 0.1733 

0.3924 


0 

0.0681 

0 



0.0487 

0.1051 



0 

0.220 



- 0.0602 

0.3220 



- 0.8029 

0.3556 



- 0.1002 

0.3808 



- 0.1064 

0.3892 



N 


1.0305 

1.0606 

1.0499 

1.0185 

0.9885 

0.882 

0.5781 

0.44 

0.1676 

- 0.1332 

- 0.2517 


0.9087 

0.8968 

0.8519 

0.7379 

0.440 

0.1442 

- 0.1713 


0.822 

0.7068 

0.440 

0.1232 

0 

- 0.0748 

- 0.1040 



















N 



Figure 20. — Correlation of Dimensionless Momentum Equation, Function Nj(n,S^) 
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with y = 1.4 used for the ratio of specific heats. The remaining parameters, the heat 
transfer factor k, the Prandtl number Pr, and the freestream Mach number. Moo, are 
input variables. 

1 The values of the correlation number at the stagnation and separation points are 

i 

( n5tag = --1064+.01725S^+.375S^ (73) 


2 

nsep= .06961 058-.0395712S^ +.06717244 5^ (74) 


and are shown in figure 21. 


Step 2 

I 

f 

The correlation number is calculated from equation (70) with the following equations 
j for the coefficients A, B, which are assumed to depend on the last available value of n 
I and S\\r. 

A(n,S^)=ai -c, n^-2d^ n^ 

^ (75) 

“ B (n,S^)-bi +2c| n + 3djn^ 

“J 

with . 


a I = .44 

bi = 5.56903 + 2.51388^ 

Cl = 3.195945 -7.08078^ 
d] =-6.358574- 13.647848^ 


* Equation (70) is solved by marching along the surface of the airfoil beginning at the 
stagnation point. The step size is the distance between the computational surface points. 
For the first and second points, n is obtained from 

i »^l = "stag ( 76 ) 

- "2"~^("l’Sw) MgVdx 1 2 2 
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From the third point on, the following numerical scheme is used. Equation (70) can be 
written as 

X 

F(x)= / f(xi)dxi 
o 


with 


F(x) = 


n M, 


B 


d Mg 

dx 




and 


M 


B-1 


f(x) = 




^ Using the trapezoidal rule 
n(Xi_|_i) = nj^j (i>2) follows. 


(77) 


^ The value of n is restricted to the range 


^stag ^ ^ ^ 1 0 *^sep 


(78) 


The factor 10 allows the numerical scheme to continue beyond the point of separation. 
This is done to avoid premature separation during the first cycles of the overall 
iteration procedure, 
f 

Step 3 

I The momentum thickness 6 is calculated from 


^ \f a„°(' '^e) d\ 

V dx 


(79) 


which can be derived from the definition of n and the relation between the momentum 
thicknesses 6 and (itp. 


78 



Pq 

0 = -2_l0 
Pe 


tr 


In order to calculate 0, the variables X, i^o and ao are obtained from 


^Sii / 

^ ^Su V 


Ksu= 198.6 °R 


( 80 ) 




') 


(81) 


The freestream stagnation value of the kinematic viscosity I'o follows from the sequence 
of calculations listed below 


T 


OO 



Tq-^Ksu (IA * 

»‘o 




1_ 

7-1 


aoo= 49.02 [ft/se^ 


^oo ^oo 


Re 




ft 


’'oo Pq 

^ Poo Poo 


The freestream stagnation value of the speed of sound is given by 

3q = 49.02 ^T^ [ft/sec] 

It should be noted that the dimension of the speed of sound, a, is ft/sec, since 


is used. 

Step 4 

The shape factor H is calculated from 

H = (^1.1138n+ 2.38411) 



r ft 

pR =4‘'.02 

secY^ 


1 + 1 




(82) 


(83) 


(84) 
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r ’ ■■ 

where the exponent a of the Prandtl number is obtained from 



Then the displacement thickness follows from 


5* = 0H 


Step 5 

The shear parameter 1 is needed to determine the skin friction coefficient. The value of 
Jf is computed from the following cubic equation in subroutine PROOT. 

Po+Pi® +P2*^+P3 8^-n = 0 

where 

2 

Po= 0.0715016-0.04559 8^ + 0.04871 

2 

p j = -0.088925 + 0.3898894 + 1.1 1892 

3 4 

+0.990225 + 1.219532 

2 - - 

P2= -1.227559-1.662158 8^-9.193 8^ 

3 4 

-13.197 8^-17.78815 8^ 

2 

P 3 = 0.73 1 2 + 4.32497 8^ + 1 2.25 1 8^ 

3 4 

+21.8919 8^+ 30.92805 8^ 

The shear parameter Jf is chosen as the lowest real value in the range 

0<K0.7 (88) 

with the assumption that 1 = 0 if no real value exists in that range. Figure 22 shows 
plots of the shear parameter for various values of Sw- 

Two different skin friction coefficients, ' Cf . and icf., are computed in subroutine 
LAMNA. They are defined by - 



80 



I("SJ 
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and 


_ _ qe 

‘^fj ^f:q 
1 1^00 


(90) 


i Here, and qoo denote the local dynamic pressure and the freestream dynamic ^ 
f pressure, respectively. The skin friction coefficient Cf. is printed, whereas Cf’; is used in 


the load calculations. 




2« 


h 

1 dMg 

i n 

Mg dx 


(91) 


qe 

Qo< 




(92) 


In the last equation Cp is the surface pressure coefficient. The Reynolds number Rex in 
equation (91) is calculated by means of 


Reg-- 


RCx Reg Q 




(93) 


Step 6 




The boundary layer thickness is found utilizing results of the Pohlhausen theory 
(ref. 14) for incompressible boundary layers. In incompressible flow (subscript i) 


6 ,- 


1 


.2 

37 ^i ^i 


315 945 9072 

Here, Aj is the Pohlhausen parameter, defined by 


du\ 

V*' dx/i 


(94) 


and related to the shape factor Hj as shown in table 3. In calculating 8 of compressible 
boundary layers, the assumption is made that equation (94) also holds in compressible 
flow, i.e,, 
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12.0 2.250 

11.0 2.253 

10.0 

9.0 "~2.273 

8.0 [J-289 

7.8 L?-293. 

7.6 L2-297^ 

7.4 ^^301 

S..9n.,ion 

Point ► 7.052 2.308 

7.0 2.309 

6.8 2.314 

6.6 2.318 


6.4 

6.2 

6.0 

5.0 

4.0 

3.0 

2.0 
1.0 

0 

- 1.0 
- 2.0 
-3.0 
-4.0 
-5.0 
- 6.0 
-7.0 
- 8.0 
-9.0 
- 10.0 
- 11.0 

Separation 

Point -12.0 

-13.0 

- 14.0 

- 15.0 




f) 


(95) 




where the compressible parameter A is obtained from the compressible shape factor H 
by means of table 3 assuming A(H) = AjlHi). The following restrictions are imposed 
on A: 

A =12 H<2.250 

A = -15 H>3.916 
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SYMBOLS OF THE LAMINAR BOUNDARY LAYER SECTION 


Theory 

Code 

Definition 

A.B 

A,B 

Coefficients of the momentum parameter N 

3o 

AO 

Freestream stagnation value of speed of sound 

at' 


Speed of sound at edge of boundary layer 

^OC 

AINF 

Speed of sound of freestream 

Cfj 

SG 

Skin friction coefficient based on 


CFI 

Skin friction coefficient based on qoo 

Cp 


Specific heat at constant pressure 

Cp 

CP 

Surface pressure coefficient 

H 

H 

Shape factor 

Htr 


Transformed shape factor 

h 


Enthalpy 

Ksu 

198.6 

Sutherland’s constant 

k 

XK 

Heat transfer factor 

( 

AL 

Shear parameter 

Me 

AME 

Local Mach number 


AMEIF 

Freestream Mach number 


(FSMCH) 


N 


Momentum parameter 

n 

CN 

Correlation number 

^sep 

SEP 

Value of n at separation 

^^stag 

CNSTAG 

Value of n at the stagnation point 



r 


Theory 

Code 

Definition 

Pr 

PRR 

Prandtl number 

P 


Static pressure 

Po 

P 

Stagnation pressure 

Pe 


Surface pressure 

Qe 


Local value of dynamic pressure 

Qoc 


Freestream value of dynamic pressure 

R 


Gas constant 

Reft 


Reynolds number per foot 


RN 

Reynolds number per foot in millions 

Rex 

REW 

Reynolds number, 

Re- 

REMOM 

Reynolds number, 

S 


Enthalpy parameter 

S\v 

SW 

Surface value of S 

T 


Static temperature 

Tw 


" . . 1 ’ rat u re 

To 

TO 

Freestream stagnation temperature in °R 

Te 


Temperature at edge of boundary layer 

Toe 

TINE 

Freestream temperature 

U 


Transformed velocity component parallel to surface 

Ue 


Value of U at the outer edge of the boundary layer 

u 


Compressible velocity component parallel to surface 

Ve 

UE 

Compressible velocity at the outer edge of the 
boundary layer 

Uoo 


Freestream velocity 

V 


Compressible velocity component noriual lo suriace 

, _ V 


Transformed velocity component normal to surface 
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Theory 

Code 

Definition 

X 

SUMS 

Arc length, surface coordinate 

X,Y 


Transformed coordinates 

x,y 


Boundary layer coordinates 

a 

ALFA 

Angle of attack in radians 

a 

ALPHA 

Exponent of Prandtl number 

y 

1.4 

Ratio of specific heats 

d 


Boundary layer thickness 


DELT 

Nondimensional boundary layer thickness 

6* 

DISP 

Displacement thickness 

8 

AMTK 

Momentum thickness 

A 

GIK 

Pohlhausen parameter 

k 

AMU 

Coefficient of viscosity temperature relation 



Dynamic viscosity 



Kinematic viscosity 


VO 

Freestream stagnation value of kinematic viscosity 

9 


Density 



Wall shear stress 


Subscripts 


e 

Outer edge of boundary layer 

stag 

Stagnation point 

i 

Incompressible 

tr 

Transformed value 

0 

Freestream stagnation value 

w 

Value at surface 

s 

Local stagnation value 

X 

Freestream 

sep 

Separation 

-4 
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TRANSITION AND LAMINAR SEPARATION (BBCAB) 


The following calculations are performed in subroutine BLTRAN, OVERLAY (4,1). 

• Natural transition location 

• Laminar separation location 

• Laminar bubble properties 

# Short bubble with turbulent reattachment 

• Long bubble with laminar separation 

• Check for user input fixed transition location 

Subroutine BLTRAN is used in conjunction with the laminar boundary layer routine 
LAMNA, which calls BLTRAN at each surface point. 

Boundary layer transition is calculated using a standard two step approach. 

1. The laminar instability point must be located, and 

2. Once the laminar boundary layer has become unstable, the search for natural 
transition is started. 

Laminar separation and stall is predicted using both Cohen-Reshotko and Curie 
criteria. Laminar bubble properties are calculated using the Goradia-Lyman criterion. 
The following paragraphs will give the details of each calculation. 

NATURAL TRANSITION 

Natural transition location is calculated using a two step approach. The first step is to 
locate the laminar instability point, i.e., that point on the surface where disturbances of 
the laminar boundary layer will be amplified. Schlichting (ref. 16) has solved the 
Orr-Sommerfeld equation, assuming Polhausen laminar velocity profiles. The results of 
this linearized stability analysis are presented in figure 23. Correlation of the 
theoretical results has been done in terms of 


versus 


Reg - 


Ue0 


dU^ 
V dx 


89 




The solid line in figure 23 is a curve fit to the stability data made by Goradia and 
given by: 

= exp ^5.46963 + 43.37458 k+ 218.28 k^- I934.6k^-23980k'^j (961 


for -.1567 « k ^.0767 

is the critical momentum thickness Reynolds number, i.e., if the local value of 
then the laminar boundary layer is stable, if 

R<=9>(Re«)crif 


the laminar boundary layer is unstable. The first point at which the laminar boundary 
layer is unstable is called the instability point. Once the instability point has been 
found, the search for the transition point begins. Granville (ref. 17) correlated 


^K)tran'K)tran-S)i 


inst 


against an average pressure gradient parameter, k, which is equal to: 


k = 


1 


s -s • 4 . 
inst s 


/ kds 


(971 


inst 


The symbols in figure 24 represent the correlation of experimental data by Granville. 
The solid line represents a curve fit to this data made by Goradia given by the 
expression, 


A 




- _2 _3 

= 825.45 + 28183.5 k + 721988 k +631 7380 k 


(98) 


'tran 

for -.05 « k « .0767. 

Thus, the first point at which 


tran 
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is called the transition point. Upon location of the natural transition point, subroutine 
LAMNA calculates the initial values of 6 and H to start the prediction of the turbulent 
flow. 


LAMINAR SEPARATION 

Laminar separation is calculated two ways: 

1. The local value of the correlation number n is compared with the separation value 
of n. If 


n > n 


sep’ 


then separation has occured according to the Cohen and Reshotko analysis. 


2. Curie writes the shear stress as: 

p i 


(99) 


where 


e2 = Fi(k)-(MU)Gi(k) 


The parameter MU is defined by 


oUe Ue 
MU =k^ - - - — 


K)' 


, dU, 

Ue~ 


U, 




dx' 


k=— U„ 


At separation 


C= o 


MU = 


Fi(k) 

G,(k) 


Thus, if the local value of MU is greater than (MU)gep > separation has occurred 
according to the Curie criteria. Goradia’s curve fits of the function Fj(k) and Gfik) 
are shown in figures 25 and 26. 

LAMINAR BUBBLE CRITERIA 

Laminar bubble criteria (long with stall or short with turbulent reattachment) are 
based on the work of (Joradia and Lyman (ref. 18). They have correlated dM^/ds with 
Reg to come up with a, critical dMg/ds to determine bubble characteristics. The following 
test is used: 





Long bubble if 


f dMe 

-02Rej-l-_>0 

j Short bubble if 

dM„ 

: -.02 Reg -I -3^ <0 

i 

! 


( 100 ) 


( 101 ) 
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SYMBOLS OF TRANSITION AND LAMINAR SEPARATION 


Theory 

Code 

Definition 

F,,G, 

Fl.Gl 

Functions of Curies’ boundary layer analysis 

H 

HMEAN 

Boundary layer shape factor 

k 

LBL 

Pressure gradient parameter 

k 

KBAR 

Average pressure gradient parameter, defined by 
equation (97) 

i 

MBL 

Parameter in Curies’ analysis 

Me 


Local Mach number 

MU 


Parameter 

n 

CN 

Correlation number of the Cohen and Reshotko 
analysis 


SEP 

Value of n at separation 

Re^ 

REMOM 

Momentum thickness Reynolds number 

S 

SUMS 

Arc length along the airfoil surface 

Ue 

UE 

Velocity at the outer edge of the boundary layer 

e 

AMTK 

Boundary layer momentum thickness 



Viscosity 

V 

XNUXl 

Kinematic viscosity 



Wall shear stress 


\ Subscripts 

f crit 

Critical value 

1 inst 

1 

Instability point 

local 

Local value 

tran 

Transition 
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TURBULENT BOUNDARY LAYER (BBCB) 


Ordinary turbulent boundary layer calculations are performed in subroutines TURBL 
and TURB of OVERLAY (4,2). 

Subroutine TURBL uses the method of Truckenbrodt with some modifications made by 
Goradia aimed at avoiding program failures in regions of flow separation. Subroutine 
TURB uses the method of Nash and Hicks, but is only employed in the last iteration 
cycle for the purpose of predicting separation. 

Method of Truckenbrodt (BBCBA) 

TruckenbrodUs turbulent boundary layer analysis is an incompressible integral method 
based on the momentum integral equation and the energy integral equation (refs. 16 
and 19 . Details of this method follow, including the modifications made by Goradia to' 
the original analysis. 


In incompressible, two-dimensional flow, the momentum integral and the energy 
integral equations read 




( 102 ) 


(103) 


pu: 


The momentum thickness ,6, and the energy dissipation thickness, 6**, are defined by 


Ug 0 = / u(Ug-u)dy 
o 

Ug u(uUu2)dy 


(104) 


(105) 


Furthermore, Ue denotes the velocity at the outer edge of the boundary layer, the 
symbol p stands for the density, and the shear stress is denoted by t having the value 
T^v at the airfoil surface. H is the ordinary shape factor, i.e., the ratio of displacement 
and momentum thicknesses. 


I 

I 


In Truckenbrodt’s analysis the shear stress integral of equation (103) is approximated 
b\ 




/ 


5 



(106) 
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and the wall shear is obtained from the Ludwieg-Tillmann formula 


w 


-.678 H/n 0\--268 


2- 

pU, 


= 0.123 10 


(¥)■ 
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Defining a second shape factor H by 


H = - 


(108) 


ar.d substituting the empirical formula equation (106) into equation (103), the energy 
ii .egral equation tak'-s the form 

u? 




(109) 


(¥)■ 


with 


z =0.0112 
1 


. ssun ing an average value of the shape factor H^v, the last equation can be integrated 
1 1 Closed form. The result is , 


e = 


n \3 + 3n 






n 1 


3 + 3n 


u: 


/ 

4 


^ 3 + 2n 


dx 


1 


1 + n 
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"l ie subscript t refers to the initial turbulent point. The coefficient A combines the 
Vi Iu< s ( f n. z. and Hav according to 


A=(l +n)^= 0.0076 
^av 


{ 111 ) 


ill addition to the quadrature formula (110) for the calculation of e the 'Pruckimbrodt 
H; tiiod utilizes a differential equation to determine the shape factor In 'I'his shape 
acioi eimation can be derived by combining the momentum integral and the energy 
ntegral (‘quations. It reads 


^ dx /U, 0\i ^ 


m 


'’w 


( 112 ) 


pU, 
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Tho two shape factors, H and H. are related by the following empirical equation 

~ 1.269 H <113) 

^ H- 0.379 

Up to this point, the described method is entirely that of Truckenbrodt. There are 
s> veral differences between equations (110), (112) and the equations programmed in 
subroutine TURBL. These differences are: 


Truckenbrodt Subroutine TURBL 

A= 0.0070 A = 0.0079 



The effects of these small changes are unknown, but they are not expected to be 
significant. 

One f)f the major areas of concern in a coupled viscous-inviscid analysis of the 
interactive type is premature boundary layer separation during the first few cycles of 
tiie iteration procedure. This problem is avoided by constraining the shape factor H to 
the range 

1.55 <H< 1.85 (114) 


\ huh. according to equation (113), corresponds to the following restrictions for the 
Siiape iactor H. 


1.21 <H<2.09 


illo) 


These constraints can also be viewed as an artificial way of modeling the flow in the 
s p.a'aied ri*gion. 


NASH AND HICKS METHOD (BBCBB) 

.he method is an integral method (ref. 20) for the prediction of ordinary turbulent 
boundary layer characteristics. It is basically incompressible, but uses simple correction 
terms to account approximately for compressibility effects. The following description 
outlines the method. 
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EQUATIONS 


':'hf Nasli and Hicks method uses the momentum integral equation 


f 9u 8u 3u , 

J J “>'1 


dy = 5 U, 


o I ' o 

a:id thi- momeni of momentum integral €>quation 


dU„ 

e dx ^ p 


1 r ^ . 

p J 9y ^ 


( I i(i) 


9u 3u 

J ^3x 3y i 

O I o 


3u , 


, 52 ., dU„ 1 p 3 t , 

y dy = Up e + — / y - 5 -“ dy 

2 e p J y dy y 
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The u-velocity is approximated by Coles’ two-parameter velocity profile iref. 8), which 
reads 


with 



K =0.41 
C- 2.05 

^ hf> two parameters in Coles' formula are the friction velocity Ux, which is related to the 
\\ .ill .-^hoar by 



ar,.: Uie velocity parameter Up 

7 nt siiear stress integral is determined from the following empirical first fjrder 
t iflareeitial equation which allows the shear stress to lag behind the equilibrium valu<‘ 
i r a L.w n value of the shape factor H. 


\vheri 



Cx 





T dy 


(119) 


1 I 2 O) 


102 



and the equilibrium value of Cjis provided by 


^ / 1 ( 121 ) 

c,= 0.025(1-^) 

hrum (•<;iiaii. ns ( 116 /, ^ 117 ^ (II81. a set of three ordinary first ord-T differemial 
euaaiio.i^ k*au he derived lor the calculation of the unknowns Uj,u^ and iht* boundary 
l.^ver tiiickness 6. The deri\ ation is outlined as follows. 


sondimensional variables a and /S are defined by 

^122) 

“=T-u; 

0 vcdocity profile, equation ( 118 ), is introduced to the momeniuni and moment 

01 momentum integral equations. This procedure provides two of the* desired ordinary 
di.feieniial ecjuations. A third equation, the so-called differentiated skin friction la'v, is 
obtained from Coles’ velocity profile by putting y = 6 in equation ( 118 ) and 
diflerentiating the result with i\‘Sj)ect to x. 


The resulting three t‘quations read 


[A]^ 


da 


dx 

d /3 


dx 

d5 


Idxj 


= {b| 


with the following coefficients of [A] and {b}. 

Aj, =-4a+ 1 - 1.58949/3 

Aj 2 = - 1.58949 a + 0.5-0.75/3 


( 123 ) 


A 13 = j/- + a - 1 .58949 oP + 0.5(J - 0.375 j 

b| =t- ^ ^4a^-3a + 3.17898aP- 1.5 p+ 0.75pV 
A2i=“-t + p(| + ;;|-0-250515) 

234[f“'-f ^“'’(■i^ 72 -«'<'™') 


- 



»2 =-d; 


^ - I ) + as ( 1 0.33402) . s - 2)(4 - 


^31 




+ C+ 1 


A 32 - 1 

*33 

Knowing tho variables a, /3, and the boundary layer momentum thickness, and the 
disj lacen.eni thickness, and the skin friction coefficient are calculated Irom 


e = 8(-2a^ + a- 1.5849 o|3 + | 
5* = 5(a + ^) 


il24) 


(1251 


Cf = 2 a 


2 2 

^ r\j^ 


(126) 


INITIAL VALUES 

T ne tui-hulent analysis is startl'd, assuming the value of the momentum tlnckm'-. o at 
h’ hisl ia;ninar point is equal to the momentum thickness at the transition point, 
n..-- a plot of will be smooth through transition, while H, c^ , and d will sho\ a 
c n in discontinuity as they go through transition. The turbulent starting proc.'ss 
o*- Sis IS oi the following steps: 

1 . The compressible laminar momentum thickness is converted to the incompressible 
value; 

»mc = «con,p(>+0-»65 ”2’' 

where is the local Mach number at transition 

*- - - — f- 
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2. The initial value of the incompressible turbulent shape factor is calculated using 
the expression 


(^initial) 


turb 


1,4754 , 

log 10 Reg 


0.9698 


(128) 


where 


Ue0 


3. 




The initial value of is calculated using the Ludwieg-Tillman formula: 

-.268 -1.561 H 


4 1 


/ \-.Z05 - 

c 0.246 fRegj e 

n.tial values of a and ^ are obtained from 


(129) 


'c, -1 


«i=V — 


and 


^1- 




- 1.58949 a H + 
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1.58949 a h)^+ 1.5 h(«(H- l)-2 ck^h) 


1 


0.75 H 


'1 V ‘qi.alion for comes from the expression for 


d*-e 


\- hi h -s solved for p. 

• . - ii.uai \aiues (actual guesses) at a and P have been obtaim-d, a Xewton 
i.a, ii.'Mn technique is used to find the values of a and p which an n o: tlie 
eq ual ions: 


a 



Reg* K 


^hc)-N!-l =0 


I 131 


O' +-^ 


h( 


2(x--a+ 1.58949 a|3-/3/2 + 0.375 /? 


2^=0 


(132) 
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A maximum of 50 iterations is permitted in the starting routine. Upon success in 
starting, new values of a, /3, Cf, and 8 are computed and the solution proceeds using the 
described predictor-corrector technique. 


^ tLl i ION METHOD 


, h p.ii'anuters a. o. and Cjare obtained from the first order ordinary diflerential 
t nations H 19' and (123i, which can be written in the form 


d0 

dx 






V .-er ■ f'* denotes any of the four unknowns. The equations are integrated using the 
following predictor-corrector technique having a maximum step size 28. The solution is 
advanced from x to x+Ax using the predictor step 


^ Ax - + 



(133) 


and the corrector step 


<^x + Ax 




COMPRESSIBILITY CORRECTION 


' 134 ) 


The Nash-Hicks analysis is incompressible, so it becomes necessary to correct the 
computed viscous flow parameters for compressibility effects. The corrections read 



( 135 ) 


(136) 


Cf = Cf. 2 

comp ‘me 1 + .065 M„ 


~ O', 


comp 


(137) 


(138) 


where Mp is the local Mach number at the edge of the boundary layer. 
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SEPARATION 


Buundarv layer separation is predicted when 

a<0.023 or /3> 1 '^39) 

Tr,,. . .quation are used in place of the theoretical values ohia.neii from 

( .il. , n 1> 01 a ^ 0 and " 1 ■ Practical experience has shown that or- .023 is mere 

reasi ahh 'n .023 conesponds to Cj < .000 18.) 
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SYMBOLS OF THE ORDINARY TURBULENT 
BOUNDARY LAYER CALCULATION * 


Theory 

Code 

Definition of Symbol 

c 

2.05 

Constant in Coles’ velocity profile 

Cf 

CF 

Skin friction coefficient 


CFI 

Shear stress integral 

A 

r - : 


C 7 


Equilibrium value of shear stress integral 

t H 

HMEAN 

Shape factor, 8*/ 6 


H 


i H 

ii\ 

Shape factor, 8**/ 6 

Me 

ML, 

Local Mach number 

f 

AME 


^ n 


Exponent 

j Reg* 


Displacement thickness Reynolds number 



Momentum thickness Reynolds number 

[ Ue 

UE 

Velocity at the outer edge of the boundary layer 



Wake component of velocity profile 

'«T 


Friction velocity 

u 


Velocity component in x -direction 

x,y 

xtempTytemp 

Boundary layer coordinates 

a 

M i 

ALPHAI 

Uy/dJeK) 

f )3 

BETAI 

^/3/Ue 


DLTA 

Boundary layer thickness 

8* 

( 

DLTAS 

Displacement thickness 

} 

i 

DELSTR 


1 

8** 


Energy dissipation thickness 

e 

TTA, 

Momentum thickness 

r ' 
1 

THETA 
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Theory 

Code 

Definition of Symbol 

r" " 



i K 

0.41 

Constant in Coles’ velocity profile 

V 

XNUXI 

Kinematic viscosity 

! ^ 
1 


Density 

1 " 


Shear stress 

1 T^v 


Wall shear stress 

i ^ 


General variable 


I 

I Subscripts 

comp Compressible 

\ 

e Outer edge of boundary layer 

^ inc Incompressible 

^ t Transition 
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CORE REGION 


'J’he con* rt‘gion is that pari ol the flow field where the wake of one airfoil component is 
separated irom the boundary layer of the next downstream airfoil by a potential core, 
(fig. 27). The following four problems are solved to determine the details of the How in 
this region. 

• Wake centerline geometry 

• Wak(‘ flow 

• Boundary layer growth 

• Knd-ot-core region 

Theoretical methods ol predicting boundary layer characteristics on the* surlace of the 
aft airfoil component are described previously in this document. The re^ador sliould note 
r.a; transition of tlie boundary layei* from laminar to turbulent flow ma\ lakie nhue in 
: :u core region. This section of th(‘ document contains a detailed accouio of tlie method 
used to trace the potential flou streamline leaving the trailing edge of the upstream 
airfoil (wake cimterline'. it further contains a brief outline of the lag entrainmfmt 
metliod of fir(*(‘n ^^hicil provides the pertinent parameters of the wake f! «\\ and 
, de. cr .bes me geometric scheme for determining the downstream boundary of the core 
' re^^ioa (eno ol core). 

WAKE CENTERLINE 

The wake centerline is that potential flow streamline which is attached to the average 
trailing edge point of an airfoil component. The problem of determining the geometry of 
the wake centerline must be solved for each component of the multielement airfoil 
during each cycle of the overall iterative solution procedure. The computer code is 
contained in subroutine WAKCL of OVERLAY (1,0) and subroutines WAKEG, WAKEJ, 
WAKES, WAKET of OVERLAY (3,0). 


STREAMLINE TRACING 

Since the potential flow problem is solved on the basis of a stream function method, 
whici among other results provides the value of the stream function for each 
stagnation streamline, it is a natural approach to also use the stream function to 
calculate the position of the wake centerline. The following assumptions are made. 

• The values and locations of all vortices yj and sources o-j are known and have a 
fixed value during the calculation of the streamline. 

^ The chord length of each segment representing the wake centerline and the total 
chord length are constant. 



Wake Centerline 


Airfoil 
trailing edge 



Confluent 
boundary layer 


Boundary layer 


Potential core 


End of the 
core region 


'Airfoil surface 


Figure 27. — Core Region 


% The source distribution simulating the displacement effect of the viscous wake 
occupies either part or the whole length of the wake centerline. 


The wake centerline is discretized as shown in figure 28, where each segment of the 
polygon is inclined to the X^-axis of the global coordinate system at an angle The 
equation of the stream function along the wake centerline can therefore be written in 
parametric form as 




where known for each airfoil component. The problem of calculating the 

centerline coordinates is nonlinear, since in equation (140) the stream function 
ie^^ends in nonlinear fashion on these coordinates. This can be seen more clearly by 
writing the value of 'k at a field pomt (X,Z) 


^(X^)=cosaZ-sinaX+ 2 2 

m= 1 j= 1 


V 


(N+N^) 


+ 2 2 

m=l j=l 


m S 
Ky Oj 


This equation is derived in the Potential Flow section which also contains the 
definitions and illustrations of all variables. In particular, it is shown there that the 
stream function influence coefficients Ky and Ky are nonlinear functions of the field 
point coordinates (X,Z). 
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Figure 28. - Discretization of Wake Centerline 

Beginning from some assumed initial position of the wake centerline, the solution of the 
nonlinear equation (140) for the desired wake centerline geometry is arrived at by 
Iteration. During a step of this iteration procedure, incremental values M of the 
segment angles are computed from 

_\dX 30 9Z 90/ J I J I j 

..ape: cript (k 1) denoLos coefficients that are known from the pr viou-^, i..., th ' k-1 ' s: 
itei hon cycle. Th e matrl.\ 

\9X 90 9Z 90 /J 

is the so-called Jacobian, whose coefficients are determined as follows. 

The derivatives of the stream function, with respect to the coordinates X and Z, are the 
potential flow velocities 


TT. = 

*^1 9Z 


-V=^ 

1 m 


(143) 


at the corner points of the wake centerline segments due to the combined effect of the 
ireestream, and all vortices yj and sources o-j of the flow field. It should be emphasized 
at this point that, during the calculation of the centerline geometry, the wake sources 
(■main fixed along the assumed initial position of the wake centerline. 
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The derivatives dX/6^ and dZIdd represent a shift of the coordinates of the i-th corner 
point caused by a change in inclination of the j-th segment of the wake centerline. 
Denoting a segment corner point by 

Pi = (Xi, Zi), 

the desired derivatives are calculated as follows. 


(144) 


^i ^ fOl 

30j l0| 


G>i) 


30; 


- Cj sin0j 
Cj cosdj 


0 <i) 


(145) 


It should be noticed that the property expressed by equation (144) leads to a triangular 
Jacobian with obvious advantages for the economy of the streamline calculation. 

The angles are calculated from equation (142) utilizing a modification of the 
familiar Newton method, that is known in the literature as the Quasi-Newton method 
(ref. 21). 


Having computed the angles A6^j, the shift of the segment corner points parallel to the 
global X,Z-coordinates are obtained from 


('^wlm 3X, 


(146) 


INITIAL POSITION 

Two cases have to be distinguished when selecting the initial position of the wake 
centerline and its total length. The calculations are performed by subroutine WAKCL of 
OVERLAY (1.0). 

Case a 

During the first cycle of the overall iteration procedure, the initial position of the wake 
|•■■n:erl^nes is chosen as shown in figure 29 for a two-component airfoil. The first wake 
• enierline is assumed to be parallel to the surface of the adjacent airfoil component. The 
distance is 
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Figure 29. — Initial Position and Length of Wake Centerlines 



Figure 30. — Selection of Segment Corner Points of Initial Wake Centerline 


^sIot'^y^TE 

where the symbols h^\ox and t^K stand for the slot height at the upstream slot exit and 
llie trailing edge thickness of the upstream airfoil, respectively. The corner points of 
wake segments are chosen, as illustrated in figure 30, by shifting the corresponding 
airfoil surface points along the surface normal. This procedure gives a wake centerline 
which extends from the slot exit to the trailing edge of the neighboring airfoil. 

The initial position and the total length of the last wake centerline are selected by 
r' .voiding the chord length of the last airfoil by 1009^. All segments representing this 
\ ake centerline are of equal length. Their total number equals the number of segments 
on the upper surface of the last airfoil. 

Case b 

B- ‘ginning with the second cycle of the overall iteration procedure, the initial position of 
the wake centerlines is the computed position for the previous cycle. 


WAKE VELOCII Y 

Wake velocity is the potential How velocity computed at points on the wake cenUTline. 
The wake centerline velocity approximates the velocities at the outer edg(‘s ol' the 
viscous wake, which are equal, since the effect of wake curvature is n(‘gUHied. The 
in\ iS( .d wake centerline velocity is the correct inner limit of the outer potential flow 
solution, and should not be confused with the viscous flow velocity along the wake 
coot* rliro'. 

COMPRESSIBILITY EFFECT 

in applying the Karman-Tsien compressibility correction, the airfoil geometry is not 
scaled from a compressible to an equivalent incompressible geometry. This is justified 
by the theoretical result that the Karman-Tsien compressibility correction does not 
distort streamlines to any significant degree during the transformation from the 
compressible to the incompressible flow domain. For the same reason, the geometry of 
the wake centerline need not be corrected for compressibility effects. The wake velocity, 
(d ccjurse, is transformed to a compressible velocity as described in the section on the 
Potential Flow solution. 


WAKE FLOW (BBCC) 


i'he n Ui‘»a is coded in subroutines WAKEI, WAKED. WADEP oi OVERLAY <4.:L. 

The properties of the wake behind each airfoil component are calculated using a version 
of the Green method, reference 7, which is based on the following assumptions,^ 

• Wake flow is two-dimensional and incompressible 
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• Wake curvature effects are neglected 

• Thin boundary layer approximations apply to wake flow 

• Secondary effects on the turbulence structure of the wake are neglected 

LAG ENTRAINMENT METHOD OF GREEN (BBCCA) 

The following description of the Green method applies to the wake flow on one side of 
the wake centerline. 


It is an integral method formulated in terms of the following three dependent variables. 
The first variable is the momentum thickness, 0, defined by 

2 

V^e= J (U^-u)udy (147) 

O 

where U\\ is the potential flow velocity at the outer edge of the wake. The coordinate y 
i- n.easured perpendicular to the wake centerline and the symbol 5 represents the 
e iro:;i the wake centerline to the outer edge of the wake. 


Tne second variable is the shape factor H, defined by 



(148) 


whert ihv displacement thickness 8"^' has the familiar definition 


Uw5 


/ (Uw - u) dy 


a49) 


The third mam variable of the prediction method is the so-called entrainment 
ctKdih jeni cp;; defined by 

The tone \ariables e. H. and Cjr are governed by the momentum integral e(iuau<c;. lae 
m:re nnient equation, and au equation for the streamvise rate of change ol ttie 
vni:a inient coefi icient. The three equations are brien> described below. Tne 
momentum integial equation reads 


dx 


= -(H + 2) 


U, 


w 


dx 


(151 
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'T'hf* <^oordinate x is tangent to the wake centerline. The entrainment equation can be 
^ .ro!;: the ctoiinilions oi the entrainment coefTicieni Cj.;. equainm (Io(l. and 

a > i..)\ Shape* paranru ter 


Hi = 


6-5 


I . y. 


1 nk 


dx 


1 ^ 

d dH, 


( 


Cg + H, (H+ 1) 


dx J 




liie :?ha,>e iuctui' ii and ii| are related b\ an empirical eijualion. which reads 

Hi = 3.15 -0.01 (H- 1)2 


' 154 ) 


Hence, 

M (H-1)2 

^^1 1.72 + 0.02 (H- 1)3 


The equation for the streamwise rate of change of the entrainment coefficient, the 
so-called lag equation, is also an empirical equation. The equation takes into acco.im 
the inliuence ol the upstream flow history on the turbulent stresses and reads 



(156) 


The various terms in equation (156) are calculated from 

x4 

(q.024-H.2ce)ce 

0.012+ 1.2 Ce 

2 

c,.= 0.024 CE+ 1.2 CE H59 
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INITIAL VALUES 


Initial values for 0 and H are provided by the boundary layer calculations at the upper- 
and lower-surface trailing edge. An initial value for the entrainment coefficient is 
assumed to be given by its equilibrium value ( 
equations ( 154), (157). (160), (161 ), (162), and (163). - 

COMPLETE WAKE SOLUTION 


The wake flow is calculated on both sides of the wake centerline solving the differential 
equations (151), (i53), and (156) in marching fashion beginning at the trailing edge of 
the upstream airfoil. At each point of the wake centerline, the wake parameters of both 
sides are calculated before the integration procedure advances to the next point. 


The main result of the wake calculation is the total displacement thickness of the wake 


6 


* 

w 




(165) 


\ hich is the sum of 8 * of the upper side (subscript u) and lower side (subscript 1) of the 
wake. 


In addition, the distance from the wake centerline to the lower edge of the wake, 
denoted by the symbol Oj. is needed to predict the end of the core region. 5 is obtained 
from 





(166) 
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At the end of the potential core region, 6j and its corresponding value at the upper side 
of the wake, 6^, are saved as initial values for the confluent boundary layer calculation. 

END OF CORE (BBCCB) 

The method is programmed in subroutine ECORE, OVERLAY (4,3 K 

The physical boundaries of the core region are shown in figure 27. The downstream 
boundary of that region is termed the end of the core. It is defined by the normal to the 
surface of the aft airfoil which passes through the point of intersection of wake and 
boundary layer edges. This definition is consistent with the aerodynamic model of the 
confluent boundary layer for which initial values must be provided along the same 
surface normal. 

The notation used in determining the end of the core is illustrated in figure 31. It is 
assumed that properties of the boundary layer beneath the potential core are known. 
The wake flow calculation proceeds in marching fashion along the wake centerline. At 
each step of the calculation it is checked whether or not the end of the core region has 
been reached. Knowing the properties of the wake at a point P\. which include th(‘ half 
width of the wake , the following calculation is performed. 

The distance P, along the surface normal is determined as described in the geometry 
section of the program for the slot height calculation. Further, the distance d along the 
the surface normal measured from the point Pj to the edge of the wake is obtained from 

As 5j 

i (167 

As COS 7 +|6j . - 5 1 . j j sin 7 

The symbol As denotes the arc length between the points P[ and Pi-i on the wake 
centerline. The angle y is formed by the normal to the wake centerline at point P, and 
the surface normal of the aft airfoil, see figure 31. 

The end of the core region has been reached if 

d ^168) 

where is the thickness of the boundary layer at the point P^. 


120 




CONFLUENT BOUNDARY LAYER (BBCD) 


FLOW MODEL 

FLOW REGIONS 

The flow downstream of the slot of a two-element airfoil configuration consists of three 
regions, shown in figure 32. These regions are termed 

• Core region 

• Confluent boundary layer 

• Ordinary turbulent boundary layer 

In general, all these flow regions can exist above the surface of the second airfoil 
component. Their existence depends on many influencing parameters involving airfoil 
geometry and flight conditions. Most regions shown in figure 32 are expected to exist if 
the relative chord length of the second airfoil is large, the gap between the two airfoils 
and the angle of attack are such that only a relatively small potential core develops, 
and, in addition, the wake of the upstream airfoil does not entirely dominate the 
spreading of the confluent boundary layer. This flow condition is often encountered on a 
wing with a leading edge flap (slat). On the other hand, for a wing with a single 
trailing edge flap and a relatively large gap, the potential core often extends beyond the 
flap trailing edge and, consequently, only the core region develops. 

The flow models of the core region and of the ordinary turbulent boundary layer are 
described previously in this document. This section describes the model of the confluent 
boundary layer, developed by Goradia (refs. 2 and 22). 

Goradia divides the confluent boundary layer into two regions - main regions I and II. 

Main region I is that flow region immediately downstream of the potential core. The 
confluent boundary layer in main region I consists of three layers, figure 33, that are 
termed 

• Wall layer 

• Jet layer 

• Wake layer 


The wall layer is the continuation of the upstream boundary layer. The jet layer and 
wake layer represent the remainder of the inner and outer part of the viscous wake of 
the upstream airfoil component, respectively. Figure 33 shows a representative velocity 
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Core region 



Figure 32. — Flow Regions Above Surface of Two-Element Airfoil 


End of potential core 

I End of wake layer 


/ 
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profile of this region. The velocity at the outer edge of the wall layer is denoted by Um- 
The symbol U\\r denotes the velocity at the common boundary of jet layer and wake 
layer. The outer edges of the three layers are denoted by the variables 85, 83, and 84. 

Main region II is the confluent boundary layer region downstream of the point where 
the wake layer disappeared. The velocity profile of this region is that of a simple wall 
jet featuring a velocity maximum only, see figure 34 . The influence of the wake of the 
upstream airfoil component is not noticeable. At the end of main region II, the jet layer 
disappears and the confluent boundary layer degenerates into an ordinary turbulent 
boundary layer. 

Having discussed the basic confiuent boundary layer model of an airfoil with two 
components, its application to the more complex flow field above a multielement airfoil 
shall now be described. Figure 35 illustrates the flow model above a high-lift airfoil 
consisting of four airfoil components, a wing with a leading edge flap, and a 
double-slotted trailing edge flap. Above the main wing and above the surface of each of 
the two trailing edge flaps, the flow field is modeled by the described basic model, which 
in general consists of a core region, main regions I and II of the confluent boundary 
layer, and an ordinary turbulent boundary layer. In other words, it is assumed that at 
each slot exit a new flow field develops, simulated by the basic flow model. This 
representation of the flow above the surface of multielement airfoils ignores the 
detailed structure of the wakes and potential cores that might still exist at the trailing 
edge of the upstream airfoil component, i.e., it is assumed that near the trailing edge of 
each airfoil component the viscous flow has always degenerated into an ordinary 
turbulent boundary layer. 


ASSUMPTIONS AND LIMITATIONS 

The following assumptions are made in the prediction of the confluent boundary layer 

characteristics. 

• The flow is two-dimensional and incompressible. 

• The effect of surface curvature is neglected. 

• The development of the various viscous layers comprising the confluent boundary 
layer is governed by the turbulent boundary layer equations for a steady mean 
flow. This and the previous assumption imply that the static pressure is constant 
in direction normal to the surface along which the confluent boundary layer 
develops. 

• The confluent boundary layer is attached to the surface over which it develops. 

• The velocity profiles of the individual viscous layers are self-similar in each region 
of the confluent boundary layer. 
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A Start of the core region 
B Start of Main Region I 
D Start of turbulent boundary layer 



The following features of the flow model of the confluent boundary layer represent 
additional limitations of its applicability. 

• The characteristics of the predicted growth of the confluent boundary layer depend 
strongly on the empirical content of the turbulent flow model. 

• The model does not account for multiple potential cores and/or multiple wakes 
which might exist above slotted flaps. 


BASIC FLOW EQUATIONS 

The governing equations of all viscous layers of the confluent boundary layer are the 
two-dimensional incompressible turbulent boundary layer equations. They read 


^ + ^ = 0 
dx 3y 


( 169 ) 


9u , 3u _ ^ /T_\ 
“ 3x 9y ® dx 9y \p) 


( 170 ) 
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where u and v are the mean values of the velocity components in x- and y-directions. 
The symbols r and p denote the shear stress and fluid density, respectively, Ue is the 
velocity at the outer edge of the confluent boundary layer. The coordinates x,y are the 
usual boundary layer coordinates, where x is measured in the direction parallel to the 
airfoil surface and y along the surface normal. 

In order to solve these equations, different initial conditions and boundary conditions 
will be specified for each layer of the confluent boundary layer in subsequent sections of 
this document. The empirical input to these equations will also be given. 

MAIN REGION I (BBCDA) 

The equations used to predict confluent boundary layer characteristics of main region I 
are listed and discussed here. The method is coded in subroutine CONF7 of 
OVERLAY (4, 4). Three viscous layers comprise the confluent boundary layer of this 
region, the wall layer, the jet layer, and the wake layer (fig. 33). The governing 
equations of these layers are coupled and must be solved simultaneously, since the wall 
layer is not separated from the jet layer by a potential core. The velocity Um at the 
outer edge of the wall layer is obtained as part of the solution. 

The governing equations of the viscous layers of this region are the turbulent boundary 
layer equations (169) and (170). The equations of each layer are solved utilizing an 
integral method and the assumption of a self-similar velocity profile. Turbulence is 
modeled by empirical relations (refs. 2 and 22) for the following quantities. 

• Growth functions for the widths of the jet layer and the wake layer 

• Shear stress terms 

• Velocity profiles and their integrals 


WALL LAYER 

The solution of the wall layer has to 

y = 0 

The subscript 5 is used for parameters of the wall layer in main region I. The growth of 
the momentum thickness 0 ^ is governed by the momentum integral equation of the wall 
layer 

d 05 2 05 1 dUjyt Ug dUe^ 

^ 5-1 % dx 2 dx ^5H5 h5_i'^ 2 (172) 

% pu^ 


satisfy the following boundary conditions 
u = v = 0 r = T^ 
u = Um r = T(8^) 


(171) 
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where the shape factor H 5 is defined by 


«5 


(173) 


and the displacement thickness 8 ^ and momentum thickness are defined by 


U 


M 


h = J ^ (uM-u)dy 


(174) 


U 


M 


^5 = / ^ u(uM-u)dy 


(175) 


The momentum integral equation (172) can be derived by integrating the x-momentum 
equation (170) across the width of the wall layer utilizing the continuity equation (169), 
the boundary conditions (171), and the assumption of a one parameter velocity profile 


1 


il_ 

% \ «5/ 


(176) 


The empirical equation 


w 


-t( 65 ) 


^ = 1.385 Y 

pUj^ 


-45.79 - 0.91 8 H 5 + 17.21 Y- 0.743 


(177) 


with 


Y=ln 


%®5 


represents the difference of the wall shear stress t\y and the shear stress at the outer 
edge of the wall layer 7 ( 85 ) in equation (172). 

The energy dissipation thickness is computed by means of 


d 5 




dx 


=-35 


** 1 dUj^ 
5 U^dx 


+ 2 6 




1 1 


2-H5 Um dx 


(178) 


-25 


^ dU 

5 2 -Hc 2 dx 


/ 


U- 


M 


j 2 3 y 
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This equation is the energy integral equation of the wall layer. 8*’^ is defined by 



( 179 ) 


and H5 is the shape factor 



(180) 


The energy integral equation (178) can be derived by first multiplying the x- 
momentum equation (170) by the velocity component u and then integrating it across 
the wall layer. The definitions of the wall layer parameters, 85, 65, H5, 85**, H5, the 
boundary conditions (eq. (171)), the continuity equation (169), and the velocity profile 
(eq. (176)) are utilized during that derivation. 

The shear stress integral in equation (178) is given by the empirical formula 


I 


5 


5 



9y 


= 0.889 


23 - 158.7 - 0.636 H 5 + 48 . 55 Y- 1 . 82 y 2 
10 Y e 


(181) 


with 


Y= 8 n 


% ^5 


The shear stress t( 85) follows from 

’■(^5)_ ’’w 

2 2 ~ 2 

P Um P P U]yj 

in which the wall shear is obtained from 

T 16 - 1 14.6 - 1.819 H 5 + 35.68 Y - 1.365 Y^ 

= 0.943 10 Y • e 

P % 


(182) 


(183) 


and the second term is given by equation (177). 

Knowing the parameters 8** and ^5 the shape factor H5 can be calculated from its 
definition, but is restricted in the code to the range 



The ordinary shape factor H 5 of the wall layer, in turn is calculated from the empirical 
equation ^ ^ . 

He = 16.133-^ + (185) 


Finally, the thickness 85 of the wall layer in main region I is calculated using 

55 = 0.00434 + 9.492 05 


The displacement thickness 8 | follows from the definition of H 5 

55 = H5 05 

Having solved all parameters of the confluent boundary layer in main region I the wall 
shear stress is calculated from the Ludwieg-Tillmann formula 


2 

P Um 


= 0.123 10 


-0.678 H 5 


hs) 


-0.268 


and not from equation (183). 

The described formulation of the wall layer problem contains the velocity at the 
outer edge of the wall layer for which no equation is given. The missing equation is 
provided by the formulation of the jet layer problem. 

JET LAYER 

The momentum integral equation of the jet layer (fig. 33) can be derived assuming a 
self-similar velocity profile f(i 7 ), defined by 


u- U^- 


) f iv) 


_ 53 -y 
^ 63-65 

f (t?) = 1 .002 - 0. 1 64 77 - 1 .967 77 ^ + 1 .338 - 0.209 


The pertinent boundary conditions of the jet layer are 


y = 55 77=1 u = Um f(77) = 0 t = 7(65) 

y = 63 77 = 0 u = f ( 77 ) = 1 7=7(63) 
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Integrating the x-momentum equation (170) from y = S 5 to y = 63 and using the 
formulation of the jet layer and wall layer problems yields the momentum integral 
equation of the jet layer. Defining nondimensional velocities by 


U 

Um = u 


M 


u 

U 


w 


(192) 


and the width of the jet layer by the symbol 


bj - 53 - 55 


the equation takes the following form 


dUj^ db; dU„, dU^ 


where the coefficients are 


Cl - 2 bj U]^ Sj ^4 - Sj 44 bj - 2 8^3 b 

/_ _ 1 hV+1 

+ 2SMsbj (Um-Uw) +2 05 


5 K - Uw) 




C2 = -Sm4K-Uw) UM-SMS(UM-Uwf + Sm3(um-Uw) 

<^3 ' b Uy - S^|3 bj (um - U,„) + 2 5^,5 bj (O^ - - Dm 

+ 2 ^M3 ^j (^M “ ^w) " ^j 


(193) 


(194) 


“ ^w) ff " ^ ^M 4 ^j 


>1 Uw “ ^ ^M5 ^j (^M ^w) 


+ ^M 3 ^j^w(^M ^w) 


V H5 + H5 _ ■ 

Um-Uw)^ 5 (Hg- l p 


C 5 - 4 I 


(%-Uw) / j 


' w 


(H5-0^ P 


+ u 


- 2 / t( 53 ) t ( 65 A 


U 


M 


M I 2 2 

\.P P 


, . 5H^ j^) 


+ %(Um"^w) Hc- 


P% 
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-u 


M 


(um-Uw)(t^) J' 


with 


Sm3= f f (T?) d7? = 0.5644 

o"' 


(195) 


%4 ~ 1 " %3 " 0-4356 
%5 4t7 = 0.433 1 

The shear stress terms contained in the coefficient ^5 are represented by the previously 
listed empirical equations of the wall layer problem. In addition, the assumption is 
made that the shear stress at the outer edge of the jet layer is 

r(53)-0.3r^ (196) 


The width of the jet layer is calculated from the empirical growth function 



= 0.17 


U 


M 


U, 


M 


-U 
+ U 


(197) 


The formulation of the jet layer problem in main region I is completed by defining the 
contributions of the jet layer to the displacement thickness and the momentum 
thickness of the confluent boundary layer. 


* 

[1 - Um + Sj^3 j 


(198) 

i [um (1 

■ Um) + Sm3 Um (um - u^) 


Sm3 ■ 

Um) (um - Uw 

)-Sm5 (um-Uw)^] 

(199) 


WAKE LAYER 

The coupling of the equations of the jet layer and the wake layer is accomplished by the 
velocity U\v at their common boundary. A momentum integral equation of the wake 
layer governing the development of U\v can be obtained by using the following 
definition of the self-similar velocity profile g(i]). 

u = Ue-(u,-U^)g(r,) 


(200) 
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(201) 


T? = 


y -^3 

^V2~h 


g ( 7 ?) = 1 .0 1 94 - 0.450 7? - 0.2029 + 0. 1 543 77 ^ - 0.024 7 ]^ 


(202) 


The symbol yj /2 denotes the half velocity point of the wake layer, i.e., the y-location 
where 

u^K + Ue) 


The boundary conditions of the wake layer problem are 


y= 63 

li 

0 

C 

II 

c 

g ( 1 ?) = 1 

y = 54 

7? = K2 

u = U^ 

g ( 1 ?) = 0 


t = t(53) 
T = 0 


(203) 


Integration of the momentum equation (170) and utilization of most of the hitherto 
introduced equations of main region I yields the desired momentum equation of the 
wake layer. Defining the width of that layer by 

K = yv2-h (204) 


the equation takes the form 


d 


. ‘Ibw . . dbl . . JUm . . 
dx 2 dx 3 dx 4 dx 5 dx 


+ ^*6 


(205) 


with 

‘^l “ ^M1 ^w “ ^ ^w ^M2 {} ~ ^w) ~ ^M3 

<>2=Smi (i-U„)-S„ 2 (l-Uw)^ 

^3 = (i “ ^vv) “ ^M3 {* ~ (^M “ 

d 4 = 2 (1 -Ci„) 9; Hj 2 ‘^M3 (' -^w) •>) + >’j (' ‘ 

<*5 = SmI (1 - Ow) bw - 2 Sm2 (1 - 0„) ^ bw + (' - ^w) 
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Hs + l 


-S 


+ 2 


13 (> - Uw) fj (Um - Uw) + 2 (' - Uw) »5 H 5 

/ - \ I 2 H 5 - 5 H 5 + II 


u 


M 


^^6 “ ■ % 




,-, \5Hs-l r(is) 


Ws-D" pUm 


pU 


M 


r, 2 1^3), r. 
P U^J 


■% 


(‘-'^w)( H5-1) I ^ J2 ly (um) "y 


and 


-K 


%1 / g (t?) dTj = 1.1 


78 


r^2 2 

%2= J g^(T/) dTj = 0.786 (206) 


The width of the wake layer is governed by the empirical growth function 




dx 


Ue + U^ 


(207) 


The contributions of the wake layer to the displacement thickness and momentum 
thickness of the confluent boundary layer are computed from 



•’w ( 

Gl 


(208) 

II 

%1 

(i-Uw)-Sm2 ( 

1 \ 

(N 

1 

(209) 


The true width of the wake layer is 

64 - 63 = 2.5 

which follows directly from the definition of K2 



(210) 


( 211 ) 


The symbol K2 is used to indicate the value of tj, see equation (201), at the outer edge of 
the wake layer 84. 
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INITIAL VALUES 

The computer code contains two different sets of initial values of the variables in main 
region I. They are chosen depending on whether or not main region I is preceded by a 
core region. 


• Main region I is preceded by a core region: The flow parameters at the end of the 
core region are saved as initial conditions for main region I. The variables 8**gand 
H 5 are recomputed using 

55* = 1.73 05-0.00005 


^ 05 

«5 = 57 


1.63 <H 5 < 1.80 


H 5 =4.411 

^ Hr ~,2 


Hi 


In addition, the velocity U\y is obtained from 


( 212 ) 


^w ^‘^^wake^^e^ (213) 

w here Uwake (^ e) denotes the potential flow wake centerline velocity at the end 
of the core, x^. 


• Main region I is entered directly: Some of the initial values are given by the flow 
conditions at the slot exit, others are simply assumed. 


Wall layer parameters 

H5 = 1.6 

The symbol the momentum thickness of the wall layer at the slot exit. 

Jet layer parameters 

UM=l-01Ue (215) 

is the boundary layer thickness at the lower surface trailing edge of the 
upstream airfoil. 

Wake layer parameters 

U^ = 0.8Ug b^ = 0.4 5p (216) 

8p is the boundary layer thickness at the upper surface trailing edge of the 
upstream airfoil. 


1.68 04 


55 = 705 


(214) 



The empirical equation (186) for the thickness of the wall layer is replaced by 


65 = 0.00596 + 1 2.88 65 + 1 5.97 65 


( 222 ) 


JET LAYER 

The wake layer has disappeared, so that 


^w 


7 ( 64 )= 7 ( 83 ) = 0 


Hence, 


- 

U = — = 1 


d U, 


w 


dx 


= 0 


and the momentum integral equation of the jet layer reduces to 


(223) 

(224) 


dUM dbj dU, 
dx ^ 2 dx ‘^^dx 


(225) 


with 


Cl - 2 bj U^, S^4 - Sj^4 bj - 2 8^3 bj - l) 


+ 2 bi 




+ 205 H 5 


H 5 -H 

(H 5-I): 


(um-i) 


^ 2 ~- ^M4 (% ■ ^) ■ ^M5 (% - 1) ^ + S^j3 - l) 


C4 




2 H; 


(2 H5 -5H5 + 1) 05 (uM-l)^- 2 SM 4 bjUM(uM-l) 


(H5-1)2 

+ 2 S^3 bj (Ojyi - 1) ^ - Ujd bj - 2 Sj^5 bj (O^ - l) 


+ ^m3 (% - 1) - 2 (Um - 1) 05 


(H5-I) 


2 




= 4U 


(ll A ^5 ^(65) 

M 1% - V /H n 2 2 “ % 2 

(H5 - 0 p Um p Um 
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DISPLACEMENT THICKNESS AND MOMENTUM THICKNESS 

The displacement thickness 8* of the confluent boundary layer in main region I is the 
sum of the contributions of the wall layer, the jet layer, and the wake layer to 8*, 

5* = 65 + 5* + 

Similarly, the momentum thickness is calculated using 

0 = 05 + 0j + 0w 


(217) 


(218) 


MAIN REGION II (BBCDB) 

The properties of the confluent boundary layer of this region are computed in 
OVERLAY (4,4), subroutine CONF8. The formulation of the problem of main region II 
is very similar to that of main region I. The differences are: 

• The formulation does not contain equations for the wake layer, since the wake 
layer has disappeared at the end of main region L 

• Some empirical coefficients used in the wall layer and jet layer formulations are 
different. 


In the remainder of this section, only those equations that are different from the 
equations of main region I are discussed. 


WALL LAYER 


The first coefficient of each of the empirical equations (177), (181), and (183) for the 
shear stress terms is different. 


’’w- '^( 55 ) 
2 
P 


-45.79 -0.918 He + 17.21 Y -0.743 
1.234 Y e 


(219) 


/ 


5 


5 



23 -158.7 -0.636 He + 48.55Y -18.2 Y^ 
1.050 10 Y e 


C220) 


w 


2 

P% 


= 0.982 


16 -114.6 -1.819 He + 35.68 Y -1.365 Y" 
10 Y e 


( 221 ) 
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SH 5-1 

H 5-1 


"(»s) 
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p Um 


-u 


M 




r 


Further, the growth function for bj has a slightly different empirical coefficient. It now 
reads 


INITIAL VALUES 


db. 

=0.185 


UM-Ue 


(226) 


Main region II is never entered directly. It is either preceded by main region I or by the 
core region. The parameters at the end of the upstream region are saved as initial 
conditions for main region II. In addition 


55 =1.7305-0.00005 


(227) 


Hr = 16.133 +.^4.54 


1.63 <H 5 < 1.80 




~2 

H 5 


The last equation is different from equation (213) of main region I. 

NUMERICAL INTEGRATION METHOD 


The mathematical problem of the confluent boundary layer is formulated in terms of a 
set of n first-order ordinary differential equations and a number of algebraic equations. 
The differential equations can be written as 





j = 1, 2 , ... n 


(228) 


for the n unknowns y^. Initial values 



at the initial x-location Xq are assumed to be known. 


The equations are integrated numerically using a modification of the Euler method. For 
this purpose, discrete points xj (i = 0 ,l, 2 ...) are chosen that coincide with the 
computational surface points that are also used as segment corner points in the 
potential flow calculation. Assuming all^variab^s are known at the point Xj.i the 
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variables (j=l,2,...n) at the next point Xf are calculated using the following 
predictor-corrector type iteration procedure. 

The predictor step reads 



j 

yi -1 


+ Ax 



2 n \ 

^i -1 ^ - yi-i) 


(229) 


with 


Ax = Xj - xj_j 


The corrector step is 


/ j j i / 1 2 \ 

[ViJ =yi_i P \ymean’ymean’ -’ymean/ 


where 


J _i 

j . ( 

j^(k-i)-| 

ymean 2 

yi-1 + \y 

i) J 


(230) 


Having computed (N-l)-values of (k=2,...n) and also (y^if,^\he variables jr’i are 

calculated by taking the average 


j 

Vi 



(231) 


The code uses N=,6. The term dUg/Ue is approximated by 


dUe Uej - Uej.i 

Ue Uci + Uei_i 


(232) 


in the integration procedure. 

MODIFIED CONFLUENT BOUNDARY LAYER METHOD 

This section contains a description of a modification of the confluent boundary layer 
model, described previously, which was developed for the purpose of predicting 
separation of the confluent boundary layer. The major modification concerns the 
velocity profile of the wall layer. The power law profile of the wall layer is replaced by 
Coles’ two parameter velocity profile, which is known, for ordinary turbulent boundary 
layers, to be a realistic representation near the point of separation. Most of the 
empirical content of Goradia’s confluent boundary layer model is retained. The 
computer code is contained in OVERLAY (4,5), subroutines CONFIl, CONFI2, 
CONFDl, CONFD2, CONFPl, CONFP2, 
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COLES’ VELOCITY PROFILE 


To predict the point of separation of the confluent boundary layer, the power law 
velocity profile used for the wall layer is replaced by Coles’ two-parameter velocity 

profile (ref. 8). With the wake function approximated by a cosine, it reads 

y Ur ^ 




(233) 


where Uj is the friction velocity, which is defined in terms of the wall shear r\\r and 
density p. 


u 


r 



(234) 


Uyjis an unknown parameter with the dimensions of velocity. The constants k and C 
have the following values 


K =0.41 


C = 2.05 


(235) 


Furthermore, in equation (233) the symbols v and 85 denote the kinematic viscosity of 
air and the thickness of the wall layer, respectively. 

Introducing the velocity profile of equation (233) to the definition of displacement 
thickness, 6^, momentum thickness, ^5, and energy dissipation thickness, 6*|, of the 
wall layer results in the following three equations. 



3**fs3Q 2 2 2 

^5 " [t6 ■ 8 ^ % "F " ^ 

+3K3(^) u^ + 6(;^) Js. 


(238) 


where 


Kj = 1.589490 
K 2 - 0.697958 
K3= 1.846111 


(239) 
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MAIN REGION I 


The confluent boundary layer method is formulated as a set of eight first order ordinary 
differential equations governing the following unknowns: 

Um ' Velocity at the outer edge of the wall layer 

^5 Momentum thickness of the wall layer 

5**5 Energy dissipation thickness of the wall layer 

65 Thickness of the wall layer 

5*5 Displacement thickness of the wall layer 

Velocity at the outer edge of the jet layer 

Friction velocity 

u^ Parameter of Coles’ velocity profile 

The equations are the 

• Momentum integral equation of the wall layer 

• Energy integral equation of the wall layer 

• Momentum integral equation of the jet layer 

• Momentum integral equation of the wake layer 

• Equation (236) differentiated with respect to x 

• Equation (237) differentiated with respect to x 

• Equation (238) differentiated with respect to x 

• Differentiated skin friction law obtained from Coles velocity profile 

The momentum and energy integral equations are very similar to the ones described 
previously. However, they are given below in their most general form, i.e., they have 
not been specialized to any particular velocity profile chosen for the wall layer. 
Furthermore, most of the empirical content of the confluent boundary layer method of 
Goradia is still contained in those equations. 



The set of eight ordinary differential equations can be written as 




05, 55 , 65, 55 , U^, Ujj). 

Details of the equations are given below in terms of the coefficients of the matrices [A] 
and {B}. 

The momentum integral equation of the wall layer reads 


where 


^iid^ +^1257 


Aii = 


205-65 + 65 
Um 


Ai o “ 1 


d Ue 65 , ^w- 




pUm 


The shear stress term is represented by the empirical equation 


Tw“T 


2 

P U]^ 


= 1.385 Y 


-45.79 -0.918 H 5 + 17.21 Y - 0.743 Y^ 


Y=6n 


%^5 


and the shape factor 




The energy integral equation of the wall layer is 


21 dx 


^23 77 " ~®2 
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where 




A2i = 


3 65 * -2 65 + 265 


U- 


M 


A 23 - 1 


= - 2 U 


dUe 55-65 ^(65) y-05 a_ /u 

2 ■*■^ 7 ,2 dy\V^ 


e dx 


U 


M 


.2/ 

P Um o P 




The shear stress integral is given by the empirical equation 


5c . 23 -158.7 -0.636 He +48.55 Y -1.82 

5 T d/n\. ^ 3 (244) 


f ^ ) dy = 0.889 10 Y 

j 2 ayVu^; 
pUm 


O „TT._ ' ■ 'M- 

Further, the value of the shear stress at the outer edge of the wall layer follows from 


7 ( 65 ) _ 7-^ ^^- 7 ( 65 ) 

2 2 2 

P Uj^ P P % 


(245) 


in which the wall shear is obtained from 


'w 

2 

P 


= 0.943 10 Y 


16 -114.6 -1.819 H 5 + 35.68 Y -1.365 Y" 


(246) 


and the second term on the right hand side of equation (245) is given by equation (242). 
The momentum integral equation of the jet layer reads 


d U 


31 d^^ +^34^ +^35d^ +A36-j^-B3 


d 6< 


d 6< 


dU 


(247) 


The coefficients are 

A 31 = bj [2 - 4 S ^|3 + 3 Sj ^3 + 2 (u^ - Sj^ 5 j 

+ (%-Uw) ( 55 -S|) 

A 34 = (Um - Uw) 

^35 = - % (% - Uw) 
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^36 




i I 2 Sj^3 - S]^3 


- 2 (um-Uw) Sms] 


Bsi = U 


,.2 r( 63 ) ^2 


3 ■ *^e ^ bj + Um 


-U 


M 


P % P Um 

- (Um - IJw) [um - 2 Um Sm 3 * Uw ^MS + ('^M ' ^w) Sms] 3F 


with 


^M3 ” 0,5644 


Sm 5 = 0.4331 


The shear stress at the outer edge of the jet layer, 7 ( 83 ), is obtained from the 
assumption 


7(63) =0.3 


(248) 


and equation (246) for the wall shear. 

The shear stress at the outer edge of the wall layer, 7 ( 85 ), is given by equation (245) in 
conjunction with equations (246) and (242). 

The growth of the jet layer is calculated from the empirical function 


d b; 


dx 


J = 0.17 


Um-Uw 


The momentum integral equation of the wake layer reads 


d U 


41 +^44dT’ +^46 ^ =64 


d 


d 5; 


dU, 


with 

*41= K-Ue) (bj-bjSM3 + S5-«s) 
*44=UMK-Ue) 

*45 = - Um (Uw - Ue) 


(249) 


(250) 
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^46 ^M1 ^ ^ ^M2 ^ (^w ” ) ^M3 

®4 “ [^^e ^M 1 ~ ^ (^w “ ^e) ^M1 ''’ ^ (^w " ^e) ^M2j ^w dx 

- Um ^ * (Uw - Ue) [(^M - Uw) Sm3 ' Um] 

pUm ■> 

- (Uw - He) [Ue Smi + (u* “ ^e) Smj]^ 

Sm] = 1.178 S„2 = 0.786 


The growth of the wake layer follows from 


dbw 

dx 


= 0.185 


U„-U 

i^Tu 


(251) 


The following three equations are obtained by differentiating equations (236), (237), and 
(239) with respect to the downstream coordinate x. 


^51 dT^-^A54^ 


ux 


(252) 


where 





(253) 


d Ua| d 0 c d 5 c d Uj d ug 

A61 d^‘^^62^ + A64^ + ^ 67 + Agg = 0 


with 



146 



(254) 


A?i 


A A 

d Ui^ d 5 c d 5 c 


d Uj. 

dx 


+ A 


d u 
78 dx 




= 0 


with 


2 : 1 c* Q 2 U_ 

“255Uj^u^-4 65 Um — 


+ 3 65 Kj u^ + 6 55 


(?)’ 


^73 - % 


3 ^5 

^74 - - % ■^ 




A 77 j^l2U^|-3Ki Uf^ u^- 12 U^i +3 K 2 


% 

+ 6 K 3 — + 18-^ 


^78 --^5 [If “l3 "f % -3 Kj Uf^j ^ 


+ 6 K 2 — + 3 K 3 


The skin friction law is obtained from equation (233) by setting u=Um at y=S^. 


Ur / 65 u,- \ 

%=— («n-^+cj + u^ 


(255) 


Differentiation with respect to x yields the last desired equation. 


81 dx 


(256) 


where 

^81 " 1 


Ag4 - - 


n 5< 


^87 


_ /Um-^13 ^ 1 \ 

V Ut k / 


^88 


= - 1 


The reader should note, that the particular order in which equations and unknowns are 
arranged seems to be most efficient for the numerical solution of the set of ordinary 
differential equations (240). The chosen arrangement yields a coefficient matrix [A] 
which approximates a triangular matrix as closely as possible. 
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MAIN REGION II 


The wake layer does not exist in this region. Hence, 





(257) 


The number of unknowns reduces to seven, which are 

^ “ (^M’ ^5’ ^5 ’ ^5’ ^5 ’ 

The governing equations of main region II are derived from those of main region I by 
eliminating the momentum integral equation of the wake layer and rewriting the 
momentum integral equation of the jet layer using equation (257). All other equations 
remain the same. In particular, the empirical coefficients in the shear-stress terms and 
the growth function for the jet-layer thickness are not changed. Writing the set of 
ordinary differential equations as 

[A 1 If) = (b( < 258 . 


the coefficients of the matrices [A] and {B} are obtained from their counterparts in 
main region I as follows. 


II 

> 

^12-^12 

Bj = B 

21 -^21 

^23 = ^23 

B2 = B 


The coefficients of the momentum integral equation of the jet layer are 


^31 " [2 % - Ug - 4 S]^3 + 3 Ug S]y|3 + 2 (\]^ - Ug^ Sjy|5j 

+ (%-Ue) ( 85 - 85 ) 

A34 = (U]y[ - Ug) 

A 35 (% - Ug) 

B3 = bj[Ug-2UMSM3 + UgSM3 + 2 (u^ - Ug) 8 ^ 5 ] ^ - U 
- (% - ^e)[ - 2 Um Sm 3 + Ug 8^43 + ^Uj^ - Ug) S^sJ 


2 ’■( 85 ) 

M — T 
P u^ 


dx 
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with 


^45 " ^55 


A54 - A^4 


A^4 - A74 


INITIAL VALUES 

It is assumed that main region II is always preceded by main region I. Therefore, only 
initial values for the latter region need be specified. The initial values for the main 
region II calculation are simply the values of the variables at the end of main region I. 

In specifying initial values for main region I, two cases are distinguished depending on 
whether or not a potential core exists at the slot exit. A potential core exists at the slot 
exit if 

^Bl(’^o) + ^Sl <hslot 
The symbols have the meaning, 

^BL^Xq) Thickness of the boundary layer on the upper surface of the downstream 

airfoil at the slot exit 

6gi Boundary layer thickness at the lower surface trailing edge of the 

upstream airfoil 

hgiot Slot height at the slot exit 


d b: Um “ 

^ = 0.17 ■ ff 

dx Uf4 + Ug 


The coefficients of the remaining equations are 


A4I - ^51 


A44 - A54 


^46 ■ ^57 


^47 “ -^58 


^ 51 ~ ^61 


^52 "^62 


^56 ■ ^67 


'57 


'68 


^ 61=^71 


^63 ■ ^73 


^66 “ ^77 


^67 ■ ^^78 


A71 - A8i 


A74 - Ag4 


^ 76 “ ^87 


'77 


= A 


88 


All other coefficients are zero. 


149 


Case a 


^BL 


("o) 


1 < ^slot 


Main region I is preceded by a core region. Denoting the location of the end of the core 
by Xq the initial values read 


% (^e) 


^5-5bl(xc) 


^5 “ ^Bl(^c) 

(’‘e) 


6c =1.73 0 


BL 


(»e) 


^5 ■ ^BL 


0.8 (^e) 


Ur - U^i-, 


' w 

I — r 

P 


-^(«n^+c) 


w 


“)3“ 

-0.678 He, ,-0.268 


pUm 


= 0.123 10 


hs ) 


(259) 








Note that Ve(xe) is the compressible surface velocity of the downstream airfoil at the 
slot exit. U^ake is the compressible wake centerline velocity. 


In addition, initial values for the thicknesses of the jet layer and wake layer are needed. 

(«u) 

h = 

wake 


bj = (5i) 


K 


wake 


(260) 


with 


K2 = 2.5 

The symbol (6u)wake denotes the distance between the wake centerline and the upper 
edge of the wake; (Sjlwake denotes the corresponding value of the lower part of the 
wake. 

Case b 


^Bl(’^o) + ^S1 ^\\ot 
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A potential core does not exist. The computation enters main region I directly at the 
slot exit. The initial values at the slot exit Xq are 



^ 0 ) 

. ^5“^BLI 

(Xo) = 

0 g 2 ^5 = 1-68 0 

^5 ^BL 

(^ 0 ) = ^S 2 

^5 =5blI 

(xo) 

= 0.8 V, 




Ur / 

65 U 7 - \ 

u^= U^j, 

J 2 


V 

in-\-+c) 

1 

i/pUM 


\ 

/ 


(261) 


(’■o) 


w 


-0.678 H< 




= 0.123 10 


(^e,5)‘ 


-0.268 


Uvi0 






M 


The thicknesses of jet layer and wake layer are initially 


bj = 5 si 


K 2 


(262) 


where fig] and 8 p denote the boundary layer thicknesses at the lower- and upper-surface 
trailing edge of the upstream airfoil, respectively. 

CONFLUENT BOUNDARY LAYER THICKNESS 

The thickness 84 of the confluent boundary layer is obtained from 


64 - 85 + bj + K2 b^ 


K2 = 2.5 


(263) 


where bj and b\v are calculated from the empirical growth functions 


db 


dx 


i = 0.17 


%-Uw 


(264) 


d b„ 


dx 


= 0.185 


Ue-Uv 

Ue + U^ 


(265) 
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and the initial values of equations (260) and (262). By definition, the wake thickness 
bvv = 0 in main region II. 


Displacement thickness and momentum thickness of the confluent boundary layer are 
calculated using 


= + 6* + C 

0 = 05 + 0j + 0^ 


with 





and the constants 


Smi = 1.178 
Sm2 “ 0.786 
Sm3 = 0.5644 
Sm5 = 0.4331 


SKIN FRICTION AND SEPARATION 

Two definitions of the skin friction coefficient Cf are used 

T , 

Cf=- 


w 


-^V 2 

T c 


Cf = 


w 


P 


(266) 


(267) 


(268) 
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The wall shear r\y is computed from the friction velocity 


Hence 


which is only printed, and 



Cf=2u^ (Uoo=l) 

which is that skin friction coefficient used in the aerodynamic load computation. 
Separation of the confluent boundary layer is assumed to take place if 


Cf < 0.001 


( 269 ) 


( 270 ) 


( 271 ) 


is predicted. 



SYMBOLS OF THE CONFLUENT BOUNDARY LAYER 

The following list of symbols does not include the coding symbols of the modified 
method of Goradia. 


Theory 

Code 

Definition 

bj 

BJAVE 

Thickness of the jet layer 

b\\- 

BWAVE 

Thickness of the wake layer 

Cf 

CFIP 

Skin friction coefficient 

Cp 

CP 

Pressure coefficient 

f(Tj) 


Jet layer velocity profile 

g(Tj) 


Wake layer velocity profile 

H 5 

HAVE 

Wall layer shape factor, 6 * 5 / 6*5 

H 5 

HVFAVE 

Wall layer shape factor, 8 ** 5 /% 

bslot 

XH 

Slot height at slot exit 

n 


Exponent of wall layer velocity profile 


REAVE 

Wall layer Reynolds number based on momentum 
thickness 

Smi 

SMI 

Integral, defined by equation (206) 

Sm 2 

SM2 

Integral, defined by equation (206) 

Sm3 

SM3 

Integral, defined by equation (195) 

Sm4 

SM4 

I“Sm3 

S\15 

SMS 

Integral, defined by equation (195) 

tXK 

TETH 

Trailing edge thickness of airfoil 

u,v 


Components of velocity in directions parallel and 
normal to the surface of the airfoil 

Uo 

UE 

Velocity at the outer edge of the confluent boundary 
layer 
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Theory 

Code 

Definition 

Um 

UMAVE 

Velocity at the outer edge of the wall layer, 
nondimensionalized by Ue 

Um 


Velocity at the outer edge of the wall layer 

Uw 

UWIAVE 

Velocity at the outer edge of the jet layer, 
nondimensionalized by Ue 

Uw 


\ ol'K iu at tht oulor I'clge ol the jet layer 

x,y 


Boundary layer coordinates, x parallel to the surface, 
y normal to it 

X<) 

XINI 

X -location of the slot exit 

^1/2 

Y2CAVE 

Half- velocity point of the wake layer, where 
u ~ ]J2 U\\7 + Ue) 

So 

DLAVE 

Wall layer thickness 

83 

D3AVE 

Ou ier edge of the jet layer 

84 

D4AVE, 

D2AVE 

Outer edge of the confluent boundary layer 

Sp 

DELE 

Boundary layer thickness at the upper-surface 
trailing edge of the upstream airfoil component 

Ssi 

DELI 

Boundary layer thickness at the lower-surface 
trailing edge of the upstream airfoil component 

Ss2 

DELI2 

Thickness of the boundary layer at the slot exit 

8* 

DLSTAV 

]\ : -iart mrn; -if 'finnuem boundary layer 

S*o 

DSTRV-L 

Displacement thickness of the wall layer 

8*j ‘ 

DSTAVE, 

DSTiUT 

Displacement thickness of the jet layer 



-- - . .. 

8*w 

dstwave; 

DSTRWK 

I hs]-lacemenl thickness of the wake layer 

S**5 

D2SAVE 

Energy dissipation thickness of the wall layer 

-..V ^ * 
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Definition 


Nondimensional coordinate, defined by equations 
(189) and (201) 


I 




THAVE 


THETA2 


CF2BUM 


Wall layer momentum thickness 

Boundary layer momentum thickness at the slot exit 


Shear stress 


Nondimensional- wall shear stress 


* t(83)/pUm' 
\ ' - 

^(85)/pUm'' 


RTD3TW 


RTD5TW 


f I TWMIT5 


Ratio of the shear stress at the outer edge of the jet 
layer and the wall shear stress 

Ratio of the shear stress at the outer edge of the wall 
layer and the wall shear stress 

Difference between the wall shear stress and the 
shear stress at the outer edge of the wall layer 


idy^ SHRINT 


Wall layer shear integral, equation (181) 


f Subscripts 


Wall layer parameters in main regions I and II 
Outer edge of the confluent boundary layer 
Jet layer 


Wake layer 
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GLOBAL AERODYNAMIC PARAMETERS (BBD) 


This section contains a description of the calculation of the aerodynamic forces and 
moments acting on multielement airfoils. 


AERODYNAMIC COEFFICIENTS 
LIFT AND MOMENT COEFFICIENTS 


The lift coefficient cj of a multielement airfoil is calculated by integrating the pressure 
and friction forces. The calculations are performed in the global axis system. The forces 
and the moment acting on a multielement airfoil are 


Component of the force in direction of the global X-axis, termed axial 
force 


Component of the force in direction of the global Z-axis, termed normal 
force 


^^ 0.0 


Pitching moment about the origin of the global axis system, positive 
nose up. 


Corresponding force and moment coefficients are defined by 
A N 


=q 


M 


°° ‘"ref 


q 


"ref 


'"m 


o, o 


o, o 

2 

‘"ref 


(272) 


where 


^ OO 2 ^oo 

IS the dynamic pressure of the uniform freestream, and Cref denotes the reference chord 
length of the multielement airfoil. 


Both, the surface pressure pg and the wall shear stress contribute to these forces 
and moment coefficients. Their contributions are calculated by discretizing the airfoil 
geometry in exactly the same way as in the Potential Flow calculation, i.e., by replacing 
the actual airfoil surface by a polygon. The corner points of the polygon (figure 36) are 
positioned on the airfoil surface and are identical with the so-called computational 
surface points defined in the geometry section. Writing these corner points in terms of 
the global airfoil coordinates (Xj,Zi) the contributions of the surface pressure to Cg, Cp, 
and Cnijj jj read 


‘^ap 


1 

‘"ref 


Nc 

S 

m=l 



(273) 
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Corner point 


Z 


G 


i 



(‘^mo,o)p 


Nc 




‘'ref 


m=l 


Nm 

2 

i=2 


cpc[’<c(Xi-Xi.l)+Z4Zi-Zi.|)' 


(275) 


with Cpj. denoting the value of the surface pressure coefficient 


Cp- 


Ps-Poo 

Qoo 


(276) 


at the midpoint of the i-th airfoil segment (X^.Z^.). The coordinates of this point are 
given hy 


Xc4(Xi + Xi_i) z^=i (Zj + Zi.,) 


(277) 


The symbols N^. and Njp are the total number of airfoil components and the number of 
surface points of the m-th airfoil component, respectively. 
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The contribution of the wall shear to Ca,Cn, and CmQ is 

Nc Nm 


' 1 


“F Cref S 


S, S Cf (Xi-Xi_,) 


( 278 ) 


Nc Nm 


cp. 


1 

F Cref 


X 1 Cf (Zi-Zn) 

m=l i=2 ' 


(279) 


Nc Njn 


Ko)p=X 1 § %[Zc(x,-XM)-Xc(Zi-Zi,,)] 

*^ref 

In these equations Cf^. is the value of the skin friction coefficient 


(280) 


Cf^ 


w 

%o 


(281) 


at the midpoint of the i-th segment. Note, that for the purpose of computing the lift 
coefficient, the sign of Cf is reversed on the lower surface of each airfoil component. 


<=f 


i<I, 


stag) 


(282) 


Istag is the index of the stagnation point of the m-th component. 

Axial-force, normal-force and pitching-moment coefficients are obtained from 


c„ = c„ + c„ 

3 ; F 


C- = + c„ 

np np 


'*"o,o (‘"mo,o)p‘^(‘'"’o,o)p (283) 


The lift coefficient cg follows from 


Cg= Cn cosa - sina 


(284) 


where a is the angle of attack. 
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DRAG COEFFICIENT 


The drag coefficient of the airfoil is calculated using the Squire and Young formula 
(ref. 9). The drag coefficient cd of each surface of each of the Nc airfoil components is 
obtained from 




s 




(285) 


where the boundary layer momentum thickness 6 , the shape factor H, and the 
compressible potential flow velocity are given by their values at the trailing edge 
point. In the case of a confluent boundary layer the chosen momentum thickness is that 
of the wall layer only, since 0 of the outer wake portion of the confluent boundary layer 
is already represented by the upstream airfoil. The total profile drag of the high-lift 
airfoil Cd is the sum of the drag coefficients of the 2Nc surfaces. 
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OUTPUT OF COMPUTER PROGRAM 


The output format of the NASA-Lockheed multielement airfoil code is described in the 
sequence in which it is printed. The definitions of the symbols of the output are 
contained in the table at the end of this section. 

CASE INPUT 

Table 4 shows an example of the format in which the user defined input data is printed. 
The sequence of the printed data agrees with the input format described in the section 
titled Processing of User Input. An exception is the printing of the airfoil surface points, 
which is done in the following sequence. 

• First airfoil component, from the leading edge to the trailing edge 

X - Coordinates of upper-surface points 
Z - Coordinates of upper-surface points 
X - Coordinates of lower-surface points 
Z - Coordinates of lower-surface points 

• Second airfoil component, from the leading edge to the trailing edge 

X - Coordinates of upper surface points 
etc. 

Note the leading edge point of each airfoil component appears twice. 

GEOMETRY 

Two different versions of the geometry of a multielement airfoil are printed out. They 
are 


• Input surface points in user coordinates 

• Computational surface points in global coordinates. 

Tables 5 and 6 show examples of both types of the geometry printout. All surface point 
coordinates are multiplied by the factor Sp/Cj-ef. The symbols Sp and c^ef denote the 
scale factor and reference chord, respectively. Surface points of the other airfoil 
components are printed in the sequence defined by the user. 

CASE OUTPUT 

The computed results of each angle of attack-Mach number case are printed in the 
sequence described below. Note the term "iteration number” is used in the printout for 
a cycle of the iteration procedure. An iteration cycle is defined in the Iteration 
Procedure section. 
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RHEING FOUR FLEP^ENT HIGH LIFT AIRFOU 
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Table 5. — Input Geometry 
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Table 6. — Computational Geometry 
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DETAILED OUTPUT FOR ITERATION NUMBER 0 


\ 

f 


The computed potential flow and viscous flow parameters are printed in the following 
sequence. Each boundary layer and wake layer summary includes a column with the 
corresponding potential flow pressures. 

First Airfoil Component 

• Upper Surface 

Laminar boundary l ayer summary 
Turbulent boundary layer summary 

• Lower Surface 

Laminar boundary layer summary 
Turbulent boundary layer summary 

• Second Airfoil Component 

• Upper Surface 

Laminar boundary layer summary 
Turbulent boundary layer summary 

• Lower Surface 

Laminar boundary layer summary 
etc. 

• First Airfoil Component 

Wake layer summary 

• Second Airfoil Component 

Confluent boundary layer summary 

Wake layer summary 

etc. 

Tables 7, 8, 9, and 10 show examples of the three boundary layer and wake layer 
summaries. Obviously, the ordinary turbulent boundary layer results are only printed if 
transition from laminar to turbulent flow has occured. The user is reminded the Nash 
and Hicks method and the modified confluent boundary layer method are not used in 
this iteration cycle. Therefore, only the results of the Truckenbrodt method for ordinary 
turbulent boundary layers and the results of the Goradia method for confluent boundary 
layers are printed. 
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Table 7. — Laminar Boundary Layer Summary 
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Table 8. — Turbulent Boundary Layer Summary 
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Table 9. — Confluent Boundary Layer Summary 
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Table 10. — Wake Summary 
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LOAD SUMMARIES FOR ITERATION NUMBERS 0 TO 3 


Overall lift, drag, and moment coefficients as well as axial and normal force coefficients 
are printed out in a load summary. The various coefficients are defined in the section 
titled Global Aerodynamic Parameters. A sample output is shown in table 11. 

t DETAILED OUTPUT FOR ITERATION NUMBER 4 

This output type is similar to the detailed printout of the results of iteration number 0. 
The difference is a summary of the ordinary turbulent boundary layer results of the 
Nash and Hicks method which follows the printed results of the Truckenbrodt method. 
An example of the turbulent boundary layer summary of the Nash and Hicks method is 
given by table 12. Results of the modified confluent boundary layer method are 
contained in table 13. 

LOAD SUMMARY FOR ITERATION 4 

The format of this summary is identical to the one of the load summaries of previous 
iteration numbers. 

SUMMARY OF SURFACE DISTRIBUTIONS OF FLOW PARAMETERS 

A table summarizing final values of surface distributions of the most important 
potential flow and viscous flow parameters is printed at the end of each data case. An 
example is shown in table 14. 


Table 11. — Load Summary 

LOAOS SU»1f*ARY SHEET 

BOEING FOUR ELEMENT HIGH LIFT AIRFOIL 

fREESTREAM MATH NUMBER « •UOOU f ANGLE OF ATTACK • A.AOOOO DEGREES 

REYNOLDS NUMBER »FR FOOT • l.COCOO MILLION # REFERENCE CHORD ■ 2.00000 FEET 

ITERATION^, _ __ total LIFT . TOTAL DRAG TOT AL._M0M£NT^ _ AXIAL FORCE NORMAL FORCE 

NUMBER COEFFICIENT COEFFICIENT ABOUT (0*0) COEFFICIENT COEFFICIENT 

0 _ ^ 018793 703936 -*1203X5 _ 1.683395 _ 

1 1.6A1357 .0197^6 -.689085 -.111176 1.639660 

2 1.6A2396 .019753 -.689A07 -.110331 1.638761 

3 1«6A2792 0191(63 ru 66 96 46 ..-♦110531 _ lf639lAl_ 
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Table 12. - Turbulent Boundary Layer Summary of Nash and Hicks Method 




u. i: K 
ar 3 3 
3 Z Z 


O f*>i 

^ ^ O Oi 

^ O f- 


fv. 
r- 
a *> 

r- IV 

o 


a c> 

i« 


fsj ^ 

m ^ 

•O OD 

«# (SJ 


•D tH 

« i-H 

CD fO 
fO « 
r<0 


^«DOfnO'«n^ tH«v*-<OO^i 

o^i 

^ivruh^o’OOvHi 

O' O OD 


TTTTTT* ***';•!***»»•* 


CD <v> O fn ^ f*^<MOmOtf>r\jO'^«nOiV'<^r^QOfvfn 
iTOr^<TO<Dr^4>«o cf' *A iO<^^^rr)r^(\jrv<^oivi 
'^’^mirimiNJfVflVIVfNJINJINjrsj^IVCMrviVISj CVr^^ 
OOOOOOOOOOOOOOOOOOOOOOi 
OOOOOOOOOOOOOOOOOOOOOO' 


Oc. C UJ O Ol 

a. oe <c a. a, 

rt z z 
z X 3 z o 

a V) z a cj 


z a o VI • 
o u. z z ^ 


CSC at z' 


cr>*vr^ <riv^ino^ #n^^Oi«mf-i,fsj^fvf\jajo 
0iT\0^v>O(\j«#(r>^r^h«aDaDaD0mfv(Dt^<cif\ 
o^A;aJfn||^<.op^€DOOf~<r^Jm«rlr^^^»^aDO<^J 
0000000000»-C»— 

oooooooooooooooooooooo 


oooo<vi^«r4-40>cniv«c>naoisiw<o-.^iri<o^r^ 
fv»*»iriu>^ o i-ci^^^floo i^rj^ <o r*» o o ^ 'O c\j 
000 000i-^^i>4i>«^*HrvNivrvfNj(Mmmmirc 
OOOOOO OOOOOOO Oo ooooooo 

oooooooooooooooooooooo 


ajivj^O'OocD^ oD^m^cDOi-ti^op^ h^cc^ 

«-C • <\» o «r O f*> O O ^ CM sT ITS iC ^ CO Q- O IV cT) o 

O O'OOOO 000^»^»H*— » I— t «— « rj CVi fV CVJ C^ 
O OOOOOOOOOOOOO oo OOOO OOOl 
O O'OOOlOOOOOOOOOOO OOOO oooi 
• 


X 

o 


LU O 

LU M 

ad ^ 

O -J 


O O «\i| 

O O O ^iO 
o o o o o 
O O O O O 
0^000^ 


t 


\ 


! 


^ «r oir^u>!'^osoo<i^ iAaj^i\j>4 

m«T^CT'OrV'— lA o <V) o <Vl ^ CD <D 00 
■c cvcAm<«rr>i-4f^i-^0<Mi'C0 •4’ O-^ tnt^^Or-i 
r- rj o O' a> r^x\ ^ t\i cv»i^ ^ o o O' o oo<~«^ 

■# ■ - - 


o o 
^ o 


• m ^ ^ ^ 


1 


I 


O' CVI 

rv oo 


Om iv <-« f'- d> o CD >o m ^ h- ivj o m ^ cd 

(«>rnicDfv»0'f<>'«v)ir>i>4r^<\iiv^f*>*HOiD «oiv; 

•«*’^'•O«-»00O^^c<I3 CJOOiv)(Dcrsro(\jmm 
<ocoi*>ii>r^^^«- 4 iv<^<Njfnh>irkO^^ >c ^ cvj 
® CD r»- ir fi'fvivirvj’^«-if-HtHOrvi<^rvt<-H,^^ivl 


I I III I 


Ml t 


II I I i I I III t I i I I 

i I 


I 


I 


^ .- 4 ;(v p> 40 o o CO «r o^ r*^ o^irc^ «s t^mcn ir o . 
'©C^'OI^O^Ot^COOOO'^r-caiO' ^CDN-rOODf\.< 
ocD («>u>>of^<viO‘r^ lrc-fnoo^~tf>m ^r*- < 

CD a. CD K lOcrs «rmmivaj«-«^^ <h o o o o O i 

aj (\i ^ a* a^i ^ l^< at a (vi ai ac aj a* a.' ai aj at at ai ^ • 


•I 


I 


li- O SL C>i 

a Ui u. ^ 

— CD a I 

►- z a c >-i 

« 3 u H- 


i 


I 


*•> '4' »T CD ^ a oai>4foflo<OF^^a'^oaja^a'^^i 
c^'a-Ka'CT'a aifnrr. ,r«Doa<.acoco>caa“( 
a a a- a a a aa aaaaa* ia<D«ai^.a-zc 

O Oii'c^'ccfT'®*^ tca^ooraooac cDfaiDOc^^ cDOc 
»M>-t»Hi-i^aaia«<n.#^ir^ir ^ < a- a* cc oo oo oc O t 


Z X u. z >- 

3 u. a ^ 

Z Z Z ia z 4 

a <« z 

•-• UJ ir < ^ u. 

^ Qd O ad « O 

« _l I- M 

Z VY o H- Ui 

O ui z a- »>H 

■« uj >- < Z O 

^ Od ID UJ Z 

vs kD at X < 


ajO'cs OoDCD (DO O'lamoO'^ tam Oi-i^( 
a*a-^a*^oaDcoa-cDO'^c^iTvcDaa^caa. 
a- o ^ ■4' >4 mcnmiaia^ 4 4* 4«r. tao^}^ 
4>^9'<-«‘4^404040‘49' 404^01 
000 '^r 4 r 4 aiami*> 44 iAca -D 4 > ^ a c 


•Aj o> CD I 

4 <« O 


172 



Table 13. — Confluent Boundary Layer Summary, Modified Method of Goradia 
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Table 14. — Summary of Surface Distributions of Flow Parameters 
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SYMBOLS OF PRINTED OUTPUT 


Output 

Theory 

Definition _ ... 

ALPHA 


Angle of attack in degrees 



Axial force coefficient 


Cd 

Drag coefficient 

CF 

Cf 

Skin friction coefficient, note warning immediately 
following this list of symbols 



Pitching-moment coefficient 



Normal-force coefficient 

CP 

Cp 

Surface-pressure coefficient 

CREF 

Cref 

Reference length 

DELTA 

A 

Angle of rotation between the coordinate system of an 
airfoil component and the reference coordinate system 
in degrees 

DELTA/C 

8/ 

Nondimensional boundary-layer thickness 

DELTAl 

85 

Outer edge of the wall layer 

DELTAS 

83 

Outer edge of the jet layer 

DELTA 4 

84/cref 

Nondimensional value of the confluent boundary 
layer thickness 

DELS/C 

8*V 

Nondimensional boundary layer displacement thK‘klu^-,^ 

DM/(S/C) 

DU/DS 

a Is ) 

d(U,./Uoc) 
a (s ^ref ) 

Derivative of local Mach number With respect to arc 
length 

Derivative of the surface velocity with respect to arc 
length 

FSMACH 

Moc 

Freestream Mach number 

H 

H 

Shape factor; in the confluent boundary layer 
summary, H is the shape factor of the wall layer 

IC 


Indices of components in the order that their data is 
stored 
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Output 

Theory 

Definition 

ICR 


Index of reference component for each component 

IM 


Index of main component 

IPP 


Index of pivot point used in placing each component 

IPPR 


Index of pivot point on reference component used in 
placing each component 

KF 

k 

Heat transfer factor 

LTRAN 


Transition option: =0 free transition; =1 fixed 
transition 

M 

Me 

Local Mach number 

N 

n 

Exponent of power-law velocity profile (reciprocal 
value) 

NA 


Number of angles of attack 

NC 

N(. 

Number of airfoil components 

NM 


Number of Mach numbers 

NPP 


Number of pivot points for each component 

NPT 


Number of input points for each component 

NSP 


Total number of computational surface points 

PR 

Pr 

Prandtl number 

RECRIT 


Momentum thickness Reynolds number at the point 
of instability 

RTRAN 


Momentum thickness Reynolds number at transition 

RN 

Reft^O^*^ 

Reynolds number per foot in millions 

SF 


Scale factor of conversion of input geometry to feet 

S/C 

s/Cref 

Nondimensional arc length 

SCRIT/C 

(s/Cref>inst 

Location of the point of instability 

STRAN/C 


Transition location 
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Output 

Theory 

Definition 

THETA/C 

^/Cref 

Nondimensional momentum thickness; in a confluent 
boundary layer summary, this is the value of the wall 
layer only 

THETAI/C 

8 1 /Cpef 

Initial value of 61 Cj-ef of fho turbulent boundary 
layer calculation 

TO 

T. 

Freestream stagnation temperature in °R 

UE/UF 

Ue'Ux 

Ratio of the velocity at the outer edge of the confluent 
boundary layer and the freestream velocity 

UL/UE 

Um'Ue 

Ratio of the velocity at the outer edge of the wall 
layer and the velocity at the outer edge of the 
confluent boundary layer 

UW/UE 

Uw/Ue 

Ratio of the velocity at the outer edge of the jet layer 
and the velocity at the outer edge of the confluent 
boundary layer 

vrvo 

Vc/Uoc 

Ratio of compressible surface velocity and freestream 
velocity 

x/c,z/c 

X/CREF, 

Z/CREF 

^G^^ref 

^G^^ref 

Nondimensional coordinates of the global axis system 

X(M.S.), 

Z(M.S.) 

XpSp/Cref 

ZpSp/Cj-ef 

Pivot point coordinates in the global axis system 
scaled by Sp/c ref 


(Xp)jSp/Cj-0f 

(Zp)iSp/ci.ef 

Pivot point coordinates in input coordinates of an 
individual airfoil component scaled by Sp/Cref 

XTRAN, 

ZTRAN 

(X|,Zi)ti.an 

Location of the transition point in input coordinates 

XL,ZL 


1 

Lower-surface point coordinates in input axis system 

XU,ZU 

(X],Zi)u 

Upper-surface point coordinates in input axis system 

YIC 

Yj/ 2 /Cref 

Half velocity point of the wake layer 


Note: Two different definitions of the skin friction coefficients CF are used. In all 

boundary layer summaries, the skin friction is referred to the local dynamic 
pressure. In the summary of the surface distributions of flow parameters, CF 
is based on the freestream value of the dynamic pressure. 
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COMPUTED RESULTS 


This section of the document summarizes evaluation results of the new version of the 
computer program. Most details of this study are described in a supplemental document 
(ref. 6). Figure 37 shows the geometry of the analyzed airfoil configurations. 
Corresponding airfoil parameters are also in reference 6, such as gap, overlap, and flap 
settings; and the investigated flight conditions including Reynolds number, Mach 
number, and angle of attack range. 

In the following test-theory comparison, three versions of the NASA/Lockheed 
multielement program are referred to: 

VERSION A 

This is the baseline version of the computer program, made operational for negative 
overlap of neighboring airfoil components. The base line version was available from the 
NASA in June 1976. 

VERSION B 

This version is described in reference 5. It differs from version A in these areas: 


• Ordinary turbulent boundary layer flow is calculated using the method of Nash 
and Hicks. 

• Profile drag is predicted by the Squire and Young formula. 

VERSION C 

This is the version described in this document. 


TEST-THEORY COMPARISONS 


BASIC GA(W)-1 AIRFOIL 

The basic GA(W)-1 airfoil was chosen to test the program capability in predicting 
performance characteristics of single airfoils. Figures 38 and 39 contain theoretical lift, 
pitching moment, and drag curves and their comparison with the experimental data of 
McGhee and Beasley (ref. 23). Both, version A and the new program version C predict 
identical lift and moment curves, which in turn agree with measured GA(W)-1 data up 
to the onset of trailing edge stall at about 8 degrees angle of attack. 
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GAl(W) - 1 






Figure 37. — Analyzed Airfoil Configurations 




Figure 38. — Lift and Pitching 
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Figure 39. — Drag Polar of GA(W)- 1 Single Airfoil 


Differences between drag polars are observed in Fig. 39. Version A, utilizing an 
integration of surface pressure and skin friction in the prediction of profile drag, gives 
the highest drag coefficients. Version C, applying the Squire and Young formula, offers 
drag values that are lower than the corresponding experimental drag coefficients. The 
lack of agreement of the three drag polars emphasizes the fact that even for single 
airfoils at low speed the problem of obtaining theoretically accurate drag data is not yet 
solved. 

GA(W)-1 WITH 30% CHORD FLAP 

The GA(W)-1 airfoil with a single 30% chord trailing edge flap served as the principal 
test case for this type of general aviation high-lift airfoil. The experimental data were 
measured by Wentz, Seetharam, and Fiscko (ref. 24 and 25). The data include global 
airfoil parameters as well as detailed surface pressures and boundary-layer 
characteristics. 

Lift- and pitching-moment characteristics of this airfoil with a flap deflection of 10 
degrees are shown in figure 40. The computed data of version C agree with the 
experimental results in the pre-stall angle of attack range, whereas, version B slightly 
mispredicts lift and moment curves. 

Differences between theoretical predictions and experimental data were noted at higher 
flap angles, but are not shown in this document. Details of these results and a 
discussion of possible reasons for the observed discrepancies are given in reference 6. 


BOEING HIGH-LIFT AIRFOIL 

The Boeing four-element high-lift airfoil, (fig. 37), was used as the main test case for 
multiple airfoils. It consists of a wing section with a leading-edge flap and a 
double-slotted trailing edge flap. Global airfoil parameters and detailed distributions of 
surface pressures and boundary layer data are available for comparisons. 

The lift and drag curves of this airfoil at a Reynolds number of two million, based on 
the wing reference chord, are given in figures 41 and 42. The experimental lift 
coefficients are balance data whereas the profile drag is obtained from wake rake 
measurements. 

All attempts failed using program version A to obtain a converged solution for this 
airfoil. Program version B arrived at converged solutions between 8 and 20 degrees 
angle of attack, but underpredicted the lift by a considerable amount, see fig. 41. The 
prediction of the lift coefficient is greatly improved by version C, but the reader should 
note that the potential flow solution already provides a very good approximation to the 
lift curve. The theoretical values of the profile drag of version C, shown in figure 42, are 
relatively close to the measured profile drag. In judging the quality of the agreement of 
the two types of drag curves, one should consider the problems of two-dimensional 
high-lift testing and the uncertainties in applying the Squire and Young formula to 
theoretical drag predictions. 




Figure 40. - Lift and Pitching Moment Characteristics of GA(W)-1 With 30% Chord Flap 








The pitching-moment characteristics of version C are compared with experimental data 
in figure 43. The discrepancy of the two curves at higher lift values is due to 
trailing-edge stall which is not modeled by the program. 

Figure 44 demonstrates the excellent convergence characteristics of the new program 
version C. 

Figures 45 and 46 contain comparisons of theoretical and experimental surface 
pressures at 8.4® angle of attack. These figures confirm the earlier findings, that 
version C indeed provides the best theoretical results. Differences between the theory of 
version C and experiment, however, are noted in the cove region of the main flap 
demonstrating the need for a model of the recirculating flow in the cove. 

Figure 47 shows boundary layer velocity profiles on the upper surface of the main 
component at several chordwise stations. The experimental velocity profiles reveal that 
very little confluence of slat wake and wing boundary layer has taken place and that an 
initially existing weak confluent boundary layer above the wing has degenerated early 
into an ordinary turbulent boundary layer. This feature of the flow field is very well 
simulated by version C. 
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CONCLUSIONS 


The aerodynamic model of the NASA-Lockheed multielement airfoil program has been 
extensively modified and most of the computer code has been rewritten using a 
structured approach to computer software design. The new version of the computer 
program has been documented in great detail and has been evaluated by comparing its 
theoretical predictions with recent experimental data of high-lift airfoils. Based on a 
relatively short evaluation phase of two months, the following conclusions about the 
reliability and quality of the program predictions are drawn. 

• The reliability of the program executions has been greatly improved. All test cases 
run have produced converged solutions within a few iteration cycles. This 
improvement is a consequence of the application of the structured approach to 
computer programming where much attention was paid to the functional 
decomposition of the aerodynamic model, its numerical implementation, and the 
data flow within the code. 

• The accuracy of the program predictions has been improved. This is due to several 
major modifications of the aerodynamic model - above all, due to the different 
representation of the viscous flow displacement effects and the improved model of 
the potential core region. 

• The computed results are consistent with the basic assumptions of the aerodynamic 
model. Best results are obtained in cases where most of the flow is attached to the 
airfoil’s surface, but the quality of the predictions gradually deteriorates with 
increasing trailing edge stall and cove separation. 

• The usefulness of the confluent boundary layer method of Goradia and its 
modification utilizing Coles’ velocity profile for the purpose of predicting the onset 
of confluent boundary layer separation has not yet been tested. Optimized 
configurations were chosen for most of the program evaluation with little 
confluence of wakes and boundary layers. 

• The performance of the program needs to be tested for configurations at off 
optimum shape design. 

The evaluation of the computer program was hampered by the shortage of reliable 
experimental high-lift data. Additional wind tunnel testing of some of the more 
important high-lift airfoil configurations would increase the confidence in their 
performance predictions. 

• Much additional theoretical work on two-dimensional high-lift airfoils needs to be 
done. A solution of the following two problems would immediately widen the range 
of applicability of the program. 


193 



The first one concerns a practical model of cove separation, which should be part of 
every computational method, in the high-lift area. A simple model of the recirculating 
flow in the cove region would improve the prediction of local effects such as pressures in 
the cove region; but, more important, might provide better initial values for the 
computations in the potential core region and the subsequent confluent boundary layer 
calculation. 

The second pr jblem is that of predicting the profile drag of multielement airfoils. The 
validity of th Squire and Young formula for the drag predictions of this type of airfoil, 
which repla^ ed the pressure and skin friction integration of the earlier versions of the 
program, i-^ questionable. Improvements in the drag prediction could be made by a 
better flow model of the wake behind a high-lift airfoil. This, in turn, requires 
improver lents in the simulation of near wakes and confluent boundary layers. 

Boeing Commercial Airplane Company 
T .0. Box 3707 

Seattle, Washington 98124 
December 1977 
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